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I. INTRODUCTION 


A. SYSTEM IDENTIFICATION BASICS 
System identification concerns the modeling of systems as sets of mathematical 

equations based upon the input and output of the system. [Ref. 1: pp. 3-6]. It allows 
a model to be developed when all that is available is the result of a process or the output 
of a system and not the process or the system itself. System identification is an impor- 
tant area of study. Solution of the modeling problem offers many alternatives for the 
continued study of the svstem. Among these are: 

e Nondestructive analysis of the system. 

e Simulation studies using the model. 

e Easv adaptation of the model to changing system environment. 


e Spectral analysis of the svstem. 


Modeling can simulate the system’s operation at a fraction of the cost of actual 
svstem operation. Complex operations not possible with the actual system for fear of 
damaging it or personal injury can be simulated. This can expose how the system will 
operate in adverse conditions not normally experienced. In speech processing, modeling 
the speech process has the potential for significantly reducing the amount of information 
necessary to store in order to reproduce the speech. 

The modeling process shown in Figure 1 on page 2 assumes the unknown svstem’s 
input and the output data are available for processing. In many cases, if the svstem’s 
input is unknown or data is not available, a white noise input can be used in its place. 
The modeling process uses the input and output data to find a set of parameters which 
closely approximate the operation of the system. The better the identification technique, 
the more closely the model follows the performance of the actual system. 

Many types of models are available. This thesis investigates a linear parametric 
model that can be described by difference equations. This tvpe of model lends itself well 
to simulation on a digital computer. The frequency characteristics of the svstem deter- 
mined from the parameters of these types of models are more accurate than what can 
be determined from classical means such as FFTs. This is because classical methods use 
windows which assume data beyond their extent is zero {Ref. 2: p. 173]. This is not a 


realistic assumption. Models in this category include the moving-average (MA) model, 
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Figure 1. System identification problem 


the autoregressive (AR) model, and the autoregressive-moving-average (ARMA) model. 
In the frequency domain, MA processes are characterized bv sharp nulls and smooth 
peaks and AR processes are characterized bv smooth nulls and sharp peaks. ARMA 
processes have sharp peaks and sharp nulls [Ref. 2: p. 173]. An advantage of the MA 
process 1s its inherent stabilitv. An advantage of the AR process is the large number of 
algorithms already available for modeling systems. An advantage of the ARMA process 
is that it uses far fewer parameters than either the MA or AR process alone to model a 
svstem. This satisfies the general requirement to reduce the complexity of the model. 

In addition to a large variety of models, there are two processing modes: block and 
sequential. 

Block processing uses a fixed length block of data in the parameter estimation 
process. It ignores data before and after the block. This is not a real time processing 
method because all data must be available before processing can start. Block processing 
generally involves inversions of data matrices whose sizes are on the order of 
(V+ 39) x(V + 12 where NV is the order of the AR process and .W ts the order of the 
MA process. 


to 


Sequential processing uses new data to update the parameter estimations. It starts 
by initializing an estimate of the inverse of the data covariance as as a diagonal matrix. 
It uses each new data point to update this matrix. Then it updates the parameter esti- 
mates using the updated inverse data covariance matrix. It is a real time method. The 
algorithm to implement the sequential processing method is generally more complex 
than the block method but less computationally intensive because the matrix inversions 
are not required. 

This thesis concerns only systems represented by discrete time data uniformly sam- 
pled at a sufficient rate to meet the Nyquist criteria. 

The work in this thesis assumes that the input data 1s a wide-sense stationary ran- 
dom sequence. Tests of the algorithms used a pseudorandom Gaussian input with a 


mean of zero and a variance of one. 


B. PROBLEM STATEMENT 

The purpose of this thesis 1s to develop algorithms for modeling systems as ARMA 
processes using the method of instrumental variables (IV) and a multichannel approach. 
Tests of the methods will be conducted to determine the accuracy of their results and the 
speed with which thev converge. 

The IV approach is a modification of the method of ordinarv least squares. This 
approach is developed first as a block processing case and then converted to a sequential 
processing case. Tests are conducted of only the sequential processing case. 

Using a multichannel scheme allows the input and output data of the unknown 
system to be processed separately. This reduces the sizes of the data matrices involved 
in the modeling process. Both block and sequential processing cases are formulated but 


only the block processing case 1s tested. 


C. OVERVIEW OF THESIS 

Chapter 2 is about ARMA modeling. It also presents a detailed derivation of the 
method of ordinary least squares because it forms the basis on which other modeling 
techniques depend. 

Chapter 3 presents a modified least-squares approach called the method of instru- 
mental variables. It is attractive due to its simplicity and good noise performance. 
Chapter 3 presents results of using this method on several second and third-order test 


systems. 


Chapter 4 presents a new multichannel approach to ARMA modeling. This ap- 
proach is presented in block and sequential processing forms. This chapter also presents 
several adaptations of the block form which improve its speed of convergence. 

Chapter 5 contains a summary of the thesis and hsts topics for further research. 

The appendix contains the programs used to test the sequential [V algorithm and 
the block multichannel iterative algorithm. Subroutines common to both programs are 
grouped together and listed at the end of the appendix. 


Il. ARMA MODELING 


A. ARMA PROCESSES 

Modeling as an autoregressive-moving-average (ARMA) process has the potential 
for achieving a close fit to the system using a reduced order over that which a moving 
average or an autoregressive model alone could achieve. ARMA modeling is concerned 
with finding a set of AR parameters and MA parameters which combined describe an 
ARMA process that approximates the characteristics of a target system. 

The general form of the ARMA model is shown in Figure 2 on page 6. The output 
at time m. y(v7), 18 a linear combination of past outputs and past and present inputs. 
The a, and 6, are constants referred to as tap weights. The a, parameters form the MA 
part of the ARMA model. The 5, parameters form the AR part. In equation form the 


output of the ARMA system is represented by the following difference equation: 


N M 
x) =— ) bin) +) aun — 3) (2.1) 
ix] i=0 


where N is the order of the AR part of the ARMA model and 7 is the order of the MA 
part of the ARMA model. This means the ARMA output at the current time depends 
on the last \ values of the ARMA output. The V 6, weighting parameters determine 
exactly how the new output depends on the past outputs. The J\/ a, weighting paramie- 
ters determine how the new output depends on the current and .\f — | past inputs. 


Equation (2.1) in vector form becomes: 
y=x'6 (2.2) 
where x is a (NV+ +1) x 1 vector of input and output data values given by: 
x=[-y(n—1) —y(n-2) 2. Wn —N) x(n) x(n — 1). x(n — AD] (2.3) 
and @isa(N+ 41+ 1)x 1 vector of the AR and MA tap weights given by: 


0=[), b, oe bx Ap Qy ae aval (2.4) 





Figure 2. Structure of the ARMA model 


For V+ L — 1 data points available for y and .f+ L data points available for u, we 
can write a block equation which gives the value of the output at progressive sampling 


times: 


b, 

Op bes) x/(n + 1), 
Vr = a) ha Oe 
by 

= : (2.5) 

; Ay 
ay 

y(n) x (ne) 
ayy 


The @ row in equation (2.5) is the value of } at time m + i based on output data available 
through time »+i—1 and input data available through n+ 7. The * row is identically 


equation (2.2) at ume n+i. In vector form equation (2.5) becomes: 
y =A0 (2.6) 
where @ is defined in equation (2.4): y, the vector of output values, is given by: 
y=(va-L4l) yn—-L42) ... yx)" (2.7) 


and X is a partitioned matrix with rows comprised of data vectors exactly hke equation 
(2.3) only shifted in time. At successive sampling times, when new data is obtained, data 
used to calculate the previous output shifts one column to the right. The new data fills 


in the left most y and wu columns. The matrix X is given by: 


—y(n — L) . —n—-ytl) un-L4+1) .. unm—-—pttl) 
—yn—-L+1)... -y(n—4 42) unm-—L4+2) .. un—-pt?) 
a , i ; a ; 2.8) 
—j(n — 1) = —y(n — N) u(7) - u(n — M) 


where 7 1s defined as N'+ L and pw is defined as + L. 


If the a, and b, are estimates of the true values of the AR and MA parameters, then 
the filter output will be an estimate of the true output. We use a hat over a variable (for 
example, }') to indicate an estimated value. Rewriting equation (2.6) using the estimated 


ARMA parameters results in: 
y= X60 (2.9) 
where 6 is defined as: 
7: Se oo CeeeeealL (2.10) 
and y is now the vector of estimated output values and is given by: 
y=(por— L + i yGie + 2) ean (2)16h 


Up until now, we have discussed estimating the output of a system given its input, 
past output, and an estimate of the parameters which describe it. If, however, we know 
the output and input of the system, based on these equations, we can use them to gen- 
erate a set of a, and b, which produces an ARMA output which 1s the best possible es- 
timate of the system output. Then the a, and b, will be optimal parameters for describing 


the operation of the unknown system as an ARMA process. 


B. METHOD OF ORDINARY LEAST-SQUARES 

In this thesis we use the method of ordinary least- squares as the means of finding 
the optimum set of ARMA parameters. It 1s a Well known modeling technique. It offers 
the advantage of being widely used in the scientific community for a variety of modeling 
problems. It has been applied successfully to a large number of modeling problems with 
good results and has been successfully applied to classes of problems for which other 
methods have failed. [Ref. 3: p. 4] 

To apply the method of ordinary least-squares to system identification we form the 
error between the actual system output and the estimated output generated by the 
ARMA model. This error is given by: 


= een vex (2.12) 


where y is the vector of the actual system outputs given by equation (2.7) and y is the 
vector of ARMA outputs given by equation (2.11). {Ref. 1: p. 176] 


We then let the sum of the squares of the errors at the instances of time the meas- 
urements of the data were taken become a measure of how well the estimates approxi- 
mate the true system outputs. This measure of performance, or cost function, is denoted 


J. It is written in equation form as: 


n+L 


J= ) dacl (2.13) 


i=n+l1 


Replacing the error in equation (2.13) with its equivalent expression from equation 
(2ey.7) results in: 
"9 


J=yTy + 0'XxX"0 — 20"Xy (2.14) 


Equation (2.14) shows that the performance measure is a function of the estimated pa- 
rameters. The criterion is to mininuze the measure of performance by taking its deriva- 
tive with respect to the parameter estimates and setting it equal to zero. Then equation 


(2.14) becomes: 


“py 


eS Oxk Oe (2.15) 


MY) 


Solving for 6, the parameters, gives us the result: 
ME OG (2.16) 


Equation (2.16) is the ordinary least-squares solution for the optimum ARMA pa- 
rameters. It provides the best possible description, in a least-squares sense, of the data 
source. The resulting parameters provide the closest fit to the actual input and output 
data of the system in the sense of least-squares errors. 

Equation (2.16) uses a block processing approach. The product of X7X must be 
formed and then inverted in order to calculate @. In addition to being computationally 
intensive, the estimate cannot be updated when new data becomes available without re- 
calculating (X7X)-'. A sequential update which does not require (X7X)~' to be recalcu- 
lated is presented in the next chapter in the context of the instrumental variable method 


of least-squares. 


il, INSTRUMENTAL VARIABLE METHOD OF SYSTEM 
IDENTIFICATION 


A. INTRODUCTION 

The instrumental variable (1V) method of svstem identification 1s a variation of the 
method of ordinary least-squares. Its attraction over ordinary least-squares is that there 
1s no bias in estimating the parameters when dealing with noise [Ref 4: p. 406]. Also, 
this method is known to yield consistent estimates and remains as easy to use as the 
method of ordinarv least-squares [Ref. 3: p. 119]. 

When an additive noise term 1s present in the observable output, }(7), the output 1s 


given by: 
y(n) = w(7) + (72) (32) 


Here w(7) represents the actual output of the system and +() represents the noise. When 
this noise has a non-zero mean, using the noise corrupted output to model the unknown 
system by the ordinary least-squares approach leads to inaccurate estimates of its pa- 
rameters. The parameters are referred to as biased estimates. [Ref: 3: p. 119, Ref. 1: pp. 
192-193, and Ref. 5: p. 704]. 

The IV method shown in Figure 3 on page 1] generates an estimate of the un- 
known system's output by processing the input data through an auxiliary model which 
closely approximates the unknown system. In our implementation of the IV method, 
the auxilarv model is an ARMA model. Its output is free of the noise affecting the 
unknown svstem. The IV method uses the auxiliary model output (estimate), w. to cal- 
culate the parameters of the unknown system. Therefore the IV parameter estimates are 
not biased like those generated by the method of ordinary least-squares. 

The IV method assumes the existence of a matrix Z composed of the auxiliary 


model’s input and output data which has the following two properties [Ref. 4: p. 406]: 


eee 
: a a 
Jim 7 27K=0 a3) 


where ¢€ is the error in fitting the parameter estimates to the data and 1s given by: 


AUXILIARY MODEL 


A A 
A(z)/B(z) 





Figure 3.  Mfodeling by the instrumental variable method 


SSN (3.4) 


and Q is a nonsingular square matrix. 

The first property means Z is orthogonal to the error. This leads to the cancellation 
of the bias term inherent in ordinary least-squares techniques [Ref. 4: p. S06]. The sec- 
ond property ensures the inverse of Z7X exists. Z 1s assumed to have the same structure 
and size as the data matrix X in equation (2.8). Its contents differ in that the noise 
corrupted svstem output 3({7) in X is replaced by the output of the auxiliary model s+(7) 


in Z. The new data matrix Z 1s given bv: 


—w(n-L) ... —w(n—-n tl) un—-Lt+1)... ua—ytl) 
—w(n-L+1) .. —wWn—nt+2) uln—-L+2) .. u(n—pw 42) 
_ , _ . . - 6.5) 
—w(n— 1) ia — w(n — N) u(r) = u(n — Mf) 


Ui 


Where y is defined as 1+ L and wis defined as \f + L. Comparing X in equation (2.8) 
and Z in equation (3.5), we note the substitution of w(x) for y(n). Thus, we are now 
using estimates of the true output 1(7) instead of noise corrupted samples y(7). 

To incorporate Z mto the parameter estimation process we begin with equation 


(2.12), which we rewrite as: 
y= XO +6 (3.6) 


This equation says that the estimates of the output, given by XO, differ from the actual 


outputs, y, by some fitting error ¢«. Multiplying equation (3.6) by Z’ yields: 
ZONE (3.7) 
Equation (3.3) ensures that Z7X can be inverted. Solving for @ results in: 
6 = (Z ome aZ, ON1Z. 'e (3.8) 


The (Z7X)-'Z7y term in equation (3.8) 1s the IV estimate of the parameters. It is written 


cus 
OO: ZY (3.9) 


The (Z7X)-'Z’s term in equation (3.8) represents a potential bias in the estimate. The 
first property of the Z matrix, given in equation (3.2), ensures this bias goes to zero, 


asvmptoucally. Applying this property, equation (3.8) can be rewritten as: 


= (eo (3.10) 
Equation (3.10) grves an unbiased estimate of the ARMA parameters. [Ref. 1: pp. 
192-9 al 

Other least-squares methods avoid the bias inherent in ordinary least-squares but 
they are more complicated than the IV method to implement [Ref. 3: p. 119]. Although 
this thesis does not attempt an analysis of the IV method in the presence of noise, any 
practical system identification technique must deal with noise. Hence, the attraction of 
and the desire to use the IV method. 

Equation (3.10) represents the block processing case. It assumes N+ L—1 output 


samples and ‘f+ L input samples are available. These samples are used to calculate an 


estimate of the parameters. Samples beyond this range are not included in the esti- 
mation process. Block processing involves multiplication of two L x (V4 3/4 1) ma- 
trices to form a third matrix. Then this third matrix must be inverted. This is a 
computationally intensive process. In what follows, we present a sequential algorithm 


to compute @,, which avoids matrix inversions. 


B. SEQUENTIAL LEAST-SQUARES ESTIMATION USING INSTRUMENTAL 
VARIABLES 

A sequential process for estimating the parameters of an unknown system requires 
fewer computations than a block process. In a manner similar to that presented in Hsia 
[Ref. 3: pp. 22-25] for the general least-squares case, the block IV estimation process 
described above can be converted into a sequential IV estimation process. Using the 
Sequential process also allows the coeflicients to be updated based on the new data that 
becomes available. 

The derivation of the sequential estimation procedure consists of two parts. The 
first part is the derivation of an equation to update the data matrix, Q(#+ 1), based 
on the previous data matrix, Q(m), and the new data: wm), y(m), and uim+1) 
Where m represents the iteration. The second part involves developing an equation for 
updating the estimate of the parameters, 6. {m+ 1), based on the previous estimate, 
6,,(7) , the previous data matrix, Q(»m), and the new data: wm), y(m), and 
ulm +1). 


Define the data matrix Q(m) to be: 
OG 1Z,.X,4 (3.11) 


where Z,, is given by equation (3.5) and X,, is given by equation (2.8). The property of 
equation (3.3) assures that Q exists. Note that Q(v7) includes output data available 
through m and input data available through m+1 . Since both Z, and X,, are 
mx (N + M) matrices, Q(7) will be a (N+ VW) x (N+ AD matrix. As the number of rows 
of Z and X increase to accommodate the increasing numbers of data points, the size of 


Q will remain the same. At the next sample time, i.e., at 77+ 1, the data matrix becomes: 
Q(m+1)=[Z) a:Xmai) (3.12) 


where the data matrices at m+ 1 are given by: 


Ne 


ile ee (3.13) 
z'(m+1) 
Xn 
X41 = La (3.14) 
x’ (m +1) 


and z7(m+1) and x"(m+ 1) are vectors which contain the most recent data values. 


Substituting equations (3.13) and (3.14) into equation (3.12) and expanding, results in: 


7 me 
Qim+ 1 =|(ZE zmeD].... (3.15) 
x/(n + 1) 
Expanding further vields: 
Q(m + 1) =[ZIX,, + zm + Ix’(m+ DT! (3.16) 


In equation (3.16) we see that two terms make up the new data matrix. The Z7X, term 
is all the data that was available through time sm. The z( + 1)x7(m-+ 1) term contains 
the new data. To perform the inversion, let A=ZIN~. B=zoa+ 1), C= i¥ana 


D=x"7(m+1). Then by the matrix inversion lemma: 
Qtm+1)=A7'—A7'B(C7' + DA7'B)'DA™! (3.17) 


Substituting the appropriate expressions for A, B, C, and D back into equation (3.17) 


vields the equation: 
Q(m + 1) =(ZAXm) — (ZnXm) Zim + 1) 


© (1 +x’ (m+ 1(Z2X,) zm+ IT" (3.18) 
© x"(m+ I(ZAXm) 


Substituting Q(m) for (Z7X,,)-' reduces equation (3.18) to: 


Q(m + 1) = Q(m) — Q(m)z(m + IE + x7(m + DQ(m)zim + DTT’ 


» x/(m + 1)Q(m) SFE 


This completes the first step of the derivation. Equation (3.19) expresses Q at time 
1+ 1 in terms of the old Q and the new data. The term in the brackets is a scalar. 
Computational intensity has been reduced because a large matrix does not have to be 
generated and its inverse does not have to be calculated. 

Continuing with the derivation, the estimate 0, for data available through m can 


be written as: 
6m) = (ZX) LAN m (3.20) 
The estimate 6, for data available through m+ 1 can be written as: 
a 0 -—l5T | 
6,40 + 1) = (Ze) Loa (3.21) 


Substituting equation (3.12) into equation (3.21) results in an expression for the estimate 


of the parameters in terms of the new data matrix and all the available data given by: 


84m + 1) = Qn + ZY nat (3.22) 
: sje 
Onei t+ 1) = Q(t VIZ zm4+1)}.... (3.23) 
ym + 1) 
6m + 1) = Q(t IEZZy,, + 20m + Lym $1) (3.24) 


Substituting for Q(m + 1) from equation (3.19) and expanding results in: 


04m + 1) = Q(m)Z Ym 
— Q(m)z(m + ILI + x(m + 1)Q(n)z(m + Ly x" (m - 1Q(n)Zy,, 
+ Q(m)z(m + yn + 1) (8275) 
— Q(m)z(m + EL + x0 + Q(n)z(m + 1)T7! 
© x/(m + 1L)Q(m)z(m + Ly(m + 1) 


Although somewhat lengthy, this equation has the desired form. To simplify it, its last 


two terms can be arranged into the form: 


Q(rn)a(m + 1) = [1+ x70 + Q(m)a(m + YI" + 1IQUn)z0n +: 1} (3.26) 
© y(n 4+ 1) 


The terms within the braces can be thought of as the result of a previous application of 
the matrix inversion lemma with A —a. B=, Ga and 


D =x7(m + 1)Q(m)z(mm + 1). Reversing the lemma results in: 
Q(m)z0m + DEL + x7 + 1)Q(n)z0n + LIT yn + 1) (3°23) 


Replacing the last two terms in equation (3.25) with this result gives us: 


6 Am + l)= Q(m)Z, y,, — Q(m)zim + LLL + x (m+ 1Q(m)z0m + ie 
© x/(m + IQ(m)Zy,, + Q(im)z(m + 1) (3.28) 
© CL4 x72 (m+ YNQ(m)z0m + IT y0n + 1) 


Factoring Q(m)z(m + 1) and [1 + x’(m + 1)Q(m)z(m + 1)]-! from the last two terms re- 


duces equation (3.28) further to: 


6, (m a) Q(m)Z/ y,, + Q(m)z(m + ICL + x7 (om + 1)Q(rn)z0m + ime 


Sele 
© [uin+ 1) - x7Un + NQUN)Z, ya] “— 


Substituting equation (3.11) into equation (3.20) and then equation (3.20) into equation 


(3.29) vields the final form for the update of the estimate of the parameters: 


6 Am + 1) =0,,Am) + Q(n)a(m + 1) 


2 - : n (3.30) 
© C14+xX (n+ TDQUn)zin + I) Gyn + 1) -— x’ (m+ 1)87f{m)] 


This is the desired result for updating the estimate of the parameters. Note that lke 
equation (3.19), the matrix inversion of equation (3.21) has been reduced to inversion 
of a scalar. Equation (3.30) describes the update of Gr + 1) in terms of the previous 
estimate of the parameters, 6 ,{m), the previous data matrix, Q(m), and the new data: 


wn), y(m), and ult 1). 


C. TESTING THE SEQUENTIAL INSTRUMENTAL VARIABLE ALGORITHM 
Equations (3.19) and (3.30) above comprise the sequential IV algorithm. Several 
tests of this algorithm were made using second and third-order filters as unknown sys- 
tems. Tests were run via computer simulation using the filters to generate the output 
data. A Gaussian random process with zero mean and unit variance was used as the 


input. The input was produced by IMSL subrouune GGNML. Graphs were created 


using DISSPLA. Table 1 on page 17 shows pole and zero locations as well as numera- 
tor and denominator parameters for the test filters. 


Table 1. TEST SYSTEMS FOR THE [V MODELING METHOD 


HE ST LOCATIONS | LOCATIONS AR MA 
Bee LER OF FOlsics Or ZEROS PARAMETERS | PARAMETERS 
le 


0.4454 50.228 | 0.4+31.27 
0.445 - 50.228 | O4- jl. 073 
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0.445 + 50.228 


Se 
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S72 


0.445 -j0.228 


=) 
tw 
tn 


0.6605 | | 0.0154 
0.6647 + j0.502 | 0.0462 
().6647-j0.502 3 | 0.0462 
| 0.0154 





Results of the tests are shown in graphical form in Figure 4 on page 18, Figure 5 
on page 19. and Figure 6 on page 20. Dashed lines indicate the true values of the pa- 
rameters. Solid lines are the IV method's estimates. 

For both second-order test cases shown, the algorithm converged quickly and 
produced accurate results. For the third-order test case, convergence took longer but 
the values were accurate. <A third-order system 1s more complex than a second-order 
system, so conceivably it would require more iterations to converge. The number of 1t- 
erations required is of the same order as the method of ordinary least-squares. 

Table 2 on page 21 contains the IV algorithm’s best estimates of the parameters and 
the number of iterations reyuired to converge to those estimates. It also shows the ab- 


solute and percent errors from the true parameters. 
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Figure 4. Second-order test case T2. (A) MA parameters. (B) AR parameters. 
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Second-order test case T2N. (A) MA parameters. (B) AR parameters. 
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Figure 6. Third-order test case T3. (A) MA parameters. (B) AR parameters. 
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tables COEFFICIENT ESTIMATES BY THE IV NIODELING METHOD. 


Hest PARAMETER ABSOMEVE PERCENT ITERATIONS 
rote tele ESme Ade EINKOKX ERROR 


T2 0,500 
-0.396 
0.888 
1.000 
-0.888 
0.244 


0.0154 0.0 
0.0466 + 0.0004 
O.0475 OU) 3 
0.0169 + 0.0015 
1.000 0.0 
-1.96 -).03 
Sse -0.040 


~~ 


O.437 -0.0204 





IV. SYSTEM IDENTIFICATION USING AN ITERATIVE 
MULTICHANNEL APPROACH 


A. INTRODUCTION 

This chapter presents an alternate system identification method, the iterative multi- 
channel approach. This approach differs from the IV method and the method of ordi- 
nary least-squares presented in the preceding chapters in that it processes the input and 
output data from the unknown system in separate channels. In its block processing 
form one advantage over the IV and ordinary least-squares methods is a reduction in the 
sizes of the data matrices. As a result, the computational complexity of the multichannel 
algorithm 1s on the order of A? + A? , where .V7 is the order of the MA part and NV 1s the 
order of the AR part. In contrast, the block IV and ordinary least-squares methods re- 


quire computations on tne order of (3+ A). 


B. PREVIOUS MULTICHANNEL METHODS 

Whittle [Ref. 6: pp. 129-130] was the first to develop a multichannel solution for the 
ARMA modeling problem. He sought to extend the recursive Durbin solution for estt- 
mating the parameters of a single variable autoregressive process to a multivariable 
autoregressive process. He discovered that to do this he would have to fit the data to 
two autoregressive processes simultaneously. One of the autoregressions would use 
present data samples to predict the value of the data one time step in the future. This 
is called forward prediction. The second autoregression would use present data samples 
to predict the value of the data at the previous time instant and is referred to as back- 
ward prediction. Sometime during this research, Whittle determined that if the input 
was derived from a MA scheme, making the process ARMA, then the solution would 
remain the same provided the correlations of the input used 1n the parameter estimation 
process had shifts greater than the MA scheme. Whittle’s use of the two separate and 
simultaneous autoregressions to model an ARMA process can be thought of as a 
multichannel modeling approach. 

Further work in the area of multichannel ARMA modeling was conducted by Perry 
and Parker [Ref. 7: pp. 509-510]. They started out with the ARMA problem formulation 
discussed in Chapter 2. Using the method of ordinary least-squares to minimize the 
mean square error, they found the solution for the estimate of the parameters to be the 


Wiener solution given by: 


a2 


Ry Rye PB] | Sy (4.1) 
R.. en, ar r ‘ 


u’y u'u yu 


In equation (4.1) R,, is a matrix of autocorrelations of the past outputs, R,.,, is a matrix 
of autocorrelations of the inputs, R,,, and R,,, are crosscorrelations of the input and 
output data, b is a vector of AR parameters, and a’ is a vector of MA parameters. In 
addition, r,, 1s a vector of autocorrelations of past output data with the current output 
and r,,, 1s a vector of crosscorrelations of input data with the current output. Bv as- 
suming the first MA parameter, a’, was known, they were able to treat it as a gain and 
factor it out of all the other MA parameters. This allowed them [Ref. 7: pp. 509-510] 


to extract the CV + 1)* row and column of equation (4.1) and rewrite the solution in the 


Sr Ry, i” = Pyy - yy | 4 2) 
R,, a a Py Py | 


In equation (4.2), a’ has been rewritten as a to indicate that the a’, term has been fac- 


form: 


tored out of all of the MA parameters. Reasoning that equation (4.2), the ARMA sol- 
ution, was a gencralization of the AR solution, they figured that it must have a recursive 
solution consisting of some combination of the Levinson-Durbin algorithm, a recursive 
solution for the AR problem. They then determined equation (4.2) was in a form similar 
to Whittle’s formulation of the problem. So thev reasoned that thev could use a form 
of Whittle’s solution to solve the ARMA modeling problem. Like Whittle, their solution 
consists of a forward and a backward autoregression. It uses two coupled lattice filters 
to process the input and output data. OfF diagonal elements of the lattice coefficient 


matrices specify the coupling points of the two lattices. 


C. ITERATIVE APPROACH TO MULTICHANNEL ARMA MODELING 

This thesis proposes another solution to the ARMA modeling problem using the 
multichannel approach. It is an iterative approach with no direct coupling of the two 
channels. However, note that there is an implicit coupling in the sense that the ARMA 
system output samples }(7) are a function of the present and past input samples u(y), 
and the past output samples. This is shown in equation (2.1). The approach proposed 
uses two separate finite-impulse-response (FIR) filters to estimate the unknown system. 
One filter estimates the AR part of the unknown system. The second estimates the MA 


part of the unknown system. Figure 7 on page 25 shows the structure of this approach. 


The derivation of the equations for this approach follows the method of ordinary least- 
squares. 

a Ut) 
B(z) 
When this signal passes through C(z) . if C(z) is close to B(z), the resulting signal X(z) 1s 





From Figure 7 on page 25, the value of the signal 1(z) is seen to be 


approximately A(z)U(z). Also from Figure 7 on page 25, the value of Z(z) 1s seen to be 
A(z)U(z) provided D(z) is close to A(z). The difference of these two signals forms the 
error Which we seek to minimize by the method of ordinarv least-squares. In minimizing 
the error we seek to drive D(z) and C(z) as close as possible to A(z) and B(z) , respec- 
uvely. 

Referring to Figure 7 on page 25, signals z(7) and x(m) are defined as the outputs 


from two FIR filters and are given by the equations: 
z(n) = dyu(a) + d,u(n — 1) + dun — 2) ++ + ayyul(n — M) = u’(n)d (4.3) 
x(n) = y(n) + qy(n — 1) Fyn — 2) t+ -- tey(n — A =" (we (4.4) 


where the vectors d and c represent the systems D(z) and C(z), respectively. The vectors 


d, u(z), c, and y(z) are given by the following equations: 


d=[d d, d .. dy)’ (4.5) 
u(7) =[u(n) ula—1) un—2) ... un— AM] (4.6) 
c=[l c, & .. cy]! (4.7) 
y(n) = Bn) x1) yr — 2) xr — ANT (4.8) 


The d parameters are estimates of the MA portion of the ARMA process. The c¢ pa- 
rameters approximate its AR portion. The vector u(7) is the input data vector of length 
M, the order of the MA part, and y() is the output data vector equal of length N, the 
order of the AR part. 


Forming the error between x and z results in: 
T G 
e(n) = z(n) — x(n) = u (n)d—y (“Je (4.9) 


The least-squares performance criterion is: 
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Figure 7. Mfultichannel modeling process 


L L 
J= > een = > E00 — x(n)]° (4.10) 
i! 


n= | 


Substituting for e(7) from equation (4.9), we can write equation (4.10) in vector form as: 


ie 


J= > =y of (4.11) 


Azl 


where we have dropped the time index, n, for convenience. Expanding equation (4.11) 


results in: 
L 


.— > aut - c'yy’c ~ 2d7uy"c (4.12) 


=| 


25 


We notice that the performance criterion is a function of both d andc. Minimizing the 
performance criterion bv differentiating with respect to the vector c and equating the 


results to zero yields: 





L Le 
= =0=0+ 2) (yh = 2) (uty (4.13) 
n=] n=!) 


Similarly, differentiating the performance criterion with respect to the vector d and 


equating the result to zero vields: 


L L 


=0=0+ >» (wut - 2» cuy7h Erie 


n=!) a 


ey 
Ney 





> 
Q. 


Solving equation (4.13) for c and equation (4.14) for d results in two equations for 


estimating the AR and MA parameters of the unknown system given by: 


c=, R,,d (4.15) 
and 

keke (4.16) 
where 


Kad) ae . . eG) 
le ont 1) YA) > ; 
R,,= uu =| ce (4.17) 
n=} 


ra) ee 


iS an estimate of the input autocorrelation matrix. The elements of R,, are computed 
using the unbiased method as follows: 


L-l 
ral) =F D_ wu - 1 iwi wee ve (4.18) 


j=0 


26 





Matrices R., R 


yy? uy? 


ucal to equation (4.17), where R,, is the estimate of the output autocorrelation matrix, 


and R,, appearing in equations (4.15) and (4.16) have structures iden- 


and R,, and R,, are estimates of the crosscorrelation matrices. Note that R,, = RZ,. The 


elements of these matrices; 7,,, 7, and r,,, are computed as follows: 
L-! 
ry) = v)vG—) for 1=0, 1, 2, 0. JN (4.19) 
EE 
ry) = chr yea l) [on — 0. ee, 28 (4.20) 
and 
L-l 
nil) = te) Oui - \mionw—O. linge. «iN (4.21) 
j=0 


Up to this point. following the standard least-squares procedure has resulted in two 
dependent or coupled equations to solve for the parameters of an unknown system 
modeled as an ARMA process. How best to solve these equations? By iteration. The 
steps of the iterative process are to 

e Calculate the correlation matrices and vectors from the available data. 
e Initiahze ec. Exploit the fact that the first parameter in ¢ 1s, c= 1. 

e Solve for d from this initial ¢ 

e Solve fore fromd. 


e Repeat beginning at the third step. 


Here is a summary of the equations in proper order for implementing the algorithm: 
¢ Compute R,, , R,, and R,, from equations (4.17) to (4.21). Note that R,, = RJ. 


e [nitialization: 


%—f1 0 ... O77 (4.22) 
Seon — UU tO-A. 

d**)) — RUIR, ce (4.23) 

6 ak, Rds (4.24) 


where X is the iteration number. 


2) 


This is a very simple algorithm. For the case where the AR and MA orders are equal. 
the correlation matrices are half the size of the block data matrices which must be gen- 
erated and inverted in the IV algorithm. 
1. Testing the Multichannel Iterative Algorithm 

We tested the algorithm by computer simulation using second and third-order 
filters as unknown systems. Table 3 shows pole and zero locations as well as MA and 
AR parameters for the test filters. Data lengths of 500 data points were used to calculate 
the autocorrelation and crosscorrelation matrices. Besides the reported cases, we tested 
the algorithm on several other second, third and nuxed-order cases. The performance 


of the algorithm in all cases that we tested was fairly uniform. 


Table 3. TEST SYSTEMS FOR ITERATIVE MULTICHANNEL MODELING 
METHOD 


TES EOC TION LOCATION AR MLA 
fel eeile is Creo ES OF ZEROS PARAMETERS | PARAMETERS 
0.445 + 30.228 0.4+ 51.273 
0.445 - 30.228 0.4 -j1.273 


0.445 + j0.228 
0.445 -j0.228 


0.6605 Q.O1354 
0.6647 + j0.5020 96 0.0462 
0.0647 -j0.5020 0 Ueto 

| 0.0] 54 





Results of the tests are shown in graphical form in Figure 8 on page 29, 
Figure 9 on page 30, and Figure 10 on page 31. Dashed lines indicate the true values 
of the parameters. Solid lines show the values the algorithm calculated. 

Table 4 on page 32 contains the algorithm’s best estimates of the parameters, 
along with the number of iterations required to converge to those estimates. It also 
shows the absolute and percent errors from the true parameters. For the second-order 
test cases, convergence to the true parameter values occurred within 20 iterations. The 


third-order test case took I4 iterations to converge to its steady-state values; however, 


these values were not the true parameter values. 
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Figure 8. Second-order test case T2. (A) MA parameters. (B) AR parameters. 
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Figure 9. Second-order test case T2N. (A) MA parameters. (B) AR parameters. 
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Figure 10. Third-order test case T3. (A) MA parameters. (B) AR parameters. 
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Table 4. PARAMETER ESTIMATES BY THE ITERATIVE MULTICHANNEL 
ALGORITHME. 
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12 0.500 
-0.393 
0.890 
1.000 
-0.889 
Rea ae 
T2s\ 1.060 Q. 20 
-0.794 0. 
0.798 0. 
1.000 0. 
-0.886 0. 
0.245 a 
T3 0.0153 -0,0001 0.6 [4 
OHS 7 + eS 5) 
0.0590 ma eS Raye 
res 7 ss a3 S6. 
0.99 -0.0] 1.0 
-1.79 + 0.20 10.035 
i267 -0.305 19.40 
-0.3086 +0,1497 32.66 


2. Stopping Parameter 


In tests of third-order svstems, the parameter estimates swung through or close 
to the true coefficient values and converged to values somewhat removed froin the true 
values. We developed a stopping parameter to flug the iteration when the estimates were 
closest to the true values. This occurs when the error is smallestaminc eminem 
Figure 7 on page 25, if D(z) 1s equal to A(z) and C(z) is equal to B(z), x and z will both 
equal A(z)U(z). The error will be zero. The further removed) and) @/=)are inom 
A(z) and B(z), respectively, the larger the error becomes. 

The stopping parameter is calculated from the difference of the z and x values 
at every iteration. After the parameter vectors d and ¢ are estimated for a particular it- 


eration, the stopping parameter is calculated from: 


ex (it) = 2y(11) — xe) = Yee — Uj dy (4.25) 


Where A is the iteration number and d.u,c and y are defined in equations (4.5) to (4.8). 
The parameter vectors d and c represent the systems D(z) and C(z), respectively. 

Figure 11 on page 34 and Figure 12 on page 35 graph the stopping parameter 
(dotted line) along with the estimated and truce values of the parameters. They show that 
when the stopping parameter is smallest, the parameters are Closest to their true values. 
Table 5 shows the improvement in the estimates of the parameters resulting from 


choosing those values when the stopping parameter is smallest. 


Table 5. PARAMETER ESTIMATES WHEN STOPPING PARAMETER IS 
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eS i PARAMETER PBSOL ETE PERCENT ITERATIONS 
— Eo tM Ie ERROR. ERROR 
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0.0478 + 0.0016 
0.0536 + 0.0074 
0.0196 + 09.0042 
0.97 -0.03 
-) $4 = (Ulli 
SES. -0.197 


-0.3672 arate 


The stopping parameter can be used in a real modeling situation because it 


comes from the data and the estimates of the parameters. Another measure of how well 
the estimates of the parameters fit the actual system is the norm of the coefficient error. 
This cannot be used in a real modeling situation however, because the values of the true 
parameters are not known. We calculated it for the test cases as a check on the appro- 
priateness of using the stopping parameter. Figure 13 on page 36 and Figure 14 on 
page 37 graph the stopping parameter (dotted line) and the norm of the coefficient error 
for test cases T2 and T3. On both graphs the two curves correspond well. Both reach 
their minimum value at the same point, the point where the estimates of the parameters 
BremclOsest to their true values. 
3. Linear-prediction of the Denominator Coefficients 
The iterative algorithm detailed in equations (4.22) to (4.24) starts by initializing 


the AR parameter estimates to: 


eNO: Ae 0) (4.26) 
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Stopping parameter example for test case T2. 


Figure 11. 
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Figure 12. Stopping parameter example for test case T3. 
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Figure 13. Stopping parameter and norm of the coefficient error for test case T2. 


where c, is Known to be | from the Z transform of the ARMA difference equation, 


equation (2.1). The Z transform is given by: 


Val tez7 toe +] = Uh t+ az taz 7 te] (4.27) 
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The initial estimates for the other AR parameters are zero. This can be far from 
their actual values. A closer estimate of the other AR parameters should result in 
quicker convergence for all parameters. A closer estimate of the AR parameters can be 
obtained by using linear-prediction techniques. Figure 15 on page 39 shows the ap- 
proach used. In Figure I5 on page 39, y(m) is the output from the unknown system. 
The system C’(z), which is represented bv the vector c’ , 1s the linear-prediction filter used 
to estimate 3(7). It uses the previous nm — .V samples of the output to generate a current 


estimate. This is given by the equation: 


36 


0.02 


0.01 


STOPPING PARAMETER 
COBRFFICIENT ERROR 





0 2 4 6 6 10 22 #14 #16 18 + «20 
ITERATIONS 


Figure 14. Stopping parameter and norm of the coefficient error for test case T3. 


$(n) = y(n = We’ (4.29) 
Where c’ is a vector of the tap weights of the autoregressive process given by: 
f tf U , iG bas Ni 
C= iegeGy, oC od (4.50) 
and y(7— 1) 1s a vector of the past output values given bv: 
ve — 1 =Ern— 1) pr —- 2). tr — NY] (4.31) 


Following least-squares techniques, we form the error between the estimate and the ac- 
tual value of the output. The sum of the squares of the errors becomes the performance 
criterion. This is differentiated with respect to the tap weights and set equal to zero. 


Solving this for the tap weights results in the equation: 


Ce =Rr (4.32) 


ee ale 


ay 


This is the standard Weiner filter solution [Ref. 8: p. 32]. It tells us that the best estimate 
of the AR parameters can be found from the correlations of the output data. The matrix 
R,, is the autocorrelation matrix of the past outputs and r,, is autocorrelation vector of 
the past outputs with the current output. In all cases tested, we did not achieve any 
significant improvement in the accuracy of the estimates of the parameters, or in the 
speed of convergence, using the straight linear prediction of equation (4.32). 

A modification to this approach, which we refer to as modified linear-prediction, 
uses correlation lags beginning on the order of the MA portion of the ARMA process. 
For example, correlations for calculating R,, for a third-order system would start at a lag 
of three and increase to a lag of five. Correlations for calculating r,, would start at a lag 
of four and increase to a lag of six. This ensures that the effect of the MA part of the 
unknown system is removed from the linear-prediction of the AR part. This modified 


method of linear-prediction is given by the equation: 


Co ry (q) TG me Ty (GP) les ag ce) 

cy ad te) ryy(q) Pe Clee Jae 2) alGio 2) 
“ (4.33) 

Chon rylgqt+p—l) nl¢+tp-2) . . ryy(Q) ryy(q + Pp) 


Where g is the order of the MA portion and p is the order of the AR portion. [Ref. 2: 
oye yd 

Figure 16 on page 40 shows the results of using modified linear-prediction with 
third-order test case T3. When comparing this graph to the estimates obtained without 
linear-prediction, shown in Figure 12 on page 35, note that the vertical axes have dif- 
ferent scales. Table 6 on page 39 lists the values of the estimates at iteration 10 and 
compares them with the true values via the absolute and percent errors. A comparison 
of Table 6 on page 39 with Table 5 on page 33, the best estimates without the use of 
modified linear-prediction, shows that modified linear prediction has significantly in- 
creased the accuracy of the AR estimates at the tenth iteration. The accuracy of the 
MA estimates remains approximately the same. The tenth iteration was chosen as the 
point to select the parameter values because in both cases this was the iteration where 


the stopping parameter had the smallest value. 
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Figure 15. Linear prediction block diagram 


Table 6. PARAMETER ESTIMATES USING NIODIFIED 
LINEAR-PREDICTION FOR INITIAL ESTIMATE OF AR PARAME- 
liek. 


ibeagy PARAMETER | ABSOLUTE Pence NT DreReAtIONS 
Peek BSN te ERROR 
0.0136 eet)? 
Q.0485 SOUS 
eet 2 + 0.003) 


0.010] (uss 
1.00 0.0 
sel UA0)! 
-0.019 
+0.0125 





D. FORMULATION OF THE SEQUENTIAL MULTICHANNEL APPROACH 
To decrease the computational intensity of updating the estimates of the AR and 


MA parameters due to new data, an algorithm to sequentially process the data has been 


oo 


PARAMETER VALUES 
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16. Parameter estimates for test case T3 using modified linear-prediction for 


initial estimate of AR parameters. 
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developed. Development begins with the performance criterion seen previously for the 
block data case in equation (4.10). From that starting point, equations are developed 
which relate new estimates of the parameters to the previous estimates and the new data. 
Separate but similar equations are developed for the MA and the AR coefficients. Due 
to the similar nature of the development of these equations, only the development of the 
equations for the AR coefficients is presented here. However, the final results for both 
the AR and MA coefficients are given. 


The performance criterion can be written as: 


J=) ei = )t - xj] = Dt —yjc]° (4.34) 
j=0 
=) 


=o 


Expanding this results in: 


n 


iA oe eee Y Anan & = 
jJ= Z} Ppemeeec V7 -C ag. Vy; Cc (4.35) 


i=Q 


where ec and y are defined in equations (4.7) and (4.8) and z 1s the scalar signal at the 
output of D. Differentiating the performance criterion with respect to ¢ and setting the 


result equal to zero yields: 


n n 


C ig P 
rae O= 0) vv a D3 (4.36) 
—o i=0 
Equation (4.36) can also be written as 
n n 
Savi= (ya (4.37) 
i=0 i=0 


Solving for the AR parameter vector results in: 
n 


n 
c= yar > an (4.38) 
i=0 


i=(0 


4] 


Because ¢ is an estimate based on data available through time n we signify this by in- 


troducing the index 7 on ¢ to vield: 


n n 


Cy = va aa (4.39) 


i i=) 


The first step in formulating the sequential algorithm is to develop an update 
equation for the estimate of the AR parameters. Applying the method presented in 
Graupe [Ref. 9: p. 124], we first define a new matrix P7! as 


n 
Py =) yar (4.40) 
i=U 


This is a matrix of the output data of the unknown system. At the previous time # — | 


this matrix is written as: 


n—-1 
— Tr 
Pi = » ya (4.41) 
uo 


By substituting equation (4.40) into equation (4.39) the estimate of the AR parameters 
can be rewritten as: 


n 


Cn = FF Z;¥; (4.42) 
i=0 


The right side of equation (4.42) needs to be converted into an expression containing the 
previous estimate of the parameters plus a correction term. It needs the past value of 


¢, which is ¢,_,. From equation (4.39) we can write: 


n—] n—] 
T: pe 
“n-1 = () yx; 2a (4.43) 
i=0 i=0 
This can be rewritten as: 
n—] n—-]) 
T. 
ar 7 0) yyy’ not (4.44) 
i=0 i=0 


Premultiplying equation (4.42) by Pz! and then separating the last term from the sum- 


mation results in: 


n—1 
a , - 
Pot, = ) 29; + Zan (4.45) 
i=0 
n=l 
Substituting for >°z,y, in equation (4.45) from its equivalent expression in equation (4.44) 
1=0 


vields: 


n—1 
Prle=() yay ena + 2nVn (4.46) 
=o 


By adding y,y/c,_, to the right side of equation (4.46) and grouping it with the summa- 
tion, the summation will range from i=Oton. In order not to affect the value on the 
right side of equation (4.46), y,yZc,_, must also be subtracted from the right-hand side, 


Which vields: 


n—-] 
= Pee & a Pee & eee 
P,,-1¢ = > yar JCn-y as Zn¥n an Yn¥nn-1 a ne n©n=1 (4.47) 
i=( 


Combining y,y/c,_, With the summation as describe above results in: 
nN 
=] wee eget eth, . 
Poic=() vis ene + Zan — Yani (4.48) 
i=0 
n 
Replacing Scy,y¥7 with its equivalent expression from equation (4.40) vields: 
i=) 


P>'c, = Pac + y,(z,— icone) (4.49) 


Premultiplving by P, results in the following equation for the update of the estimate of 


the AR parameters: 
GR eS (4.50) 


This is the desired result. It relates the estimate of the parameters at time JN to the es- 


timate at the previous time, NV — 1}, plus the new output data vector, y,, and the error at 
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time A. Error is represented bv the z, — yJc,_, term. Ihe corresponding equationmcmame 


n-} 


MA parameters 1s: 
dh dha; Oey Greet cine (4.51) 
In equation (4.51), Q is a matrix of the input data of the unknown system given by: 
n 
Q;' =) ua? (4.52) 


i=0 


Finally, we need a sequential update for P, and Q, to complete the sequential algo- 
rithm. This 1s accomplished by using a form of the matrix inversion lemma. 
By pulling the last term out of the summation, the definition of Pz! given in equation 


(4.40) can be rewritten as: 
n-] 
Pi = » vai + 9nYn (4.53) 
i=0 
Substituting for Lys its equivalent expression from equation (4.41) results in 
Pee Peep th Yala (4.54) 
Inverting both sides of the equation results 1n: 
(OP oy Ae Seal (4.55) 
LetA=P7,, B=y,,C=1,andD=y,/. Then, by the matrix inversion lemma: 
P,=A '—A'B(C| +DA'B) DA | (4.56) 
Substituting the appropriate expressions into equation (4.56) results in: 
ets) Neo) nestle) sel yeaa (4.57) 
This reduces to: 


area reer) | renee A ee (4.58) 


Using this same procedure, the update for Q, 1s: 


ad 


Q,, = Oe _ Ora oF Oe ie u.Q,_, (4.59) 


A reduction in the computational intensity has been achieved by reducing the matrix 
inversion of equation (4.55) to inversion of a scalar in equation (4.57). Inversion of 
these scalars 1s much simpler than inversion of the R,, and R,, matrices of the block 
processing case. 

The sequential multichannel algorithm is summarized below: 


e The parameter update equations: 


Cn = ny + PaYe(Zn — Yana) (4.60) 
eed Oye dean) (4.61) 
e The update equations for the P and Q matrices: 
Beep ar 43, 1 lice yell, Py (4.62) 
Qn = Qray — Qra tal + Uy; QyapYnd Un Qn (4.63) 


The reduction in computational intensity comes with a trade-off. Now the algo- 
rithm) is more complex to use. Updates are required for P and Q as well as c and d where 
before, in the block multichannel algorithm, updates were only required force andd. But, 
as in the sequential IV case, an added advantage of the sequential multichannel! algo- 


rithm is it allows updates of the estimates of the parameters based upon new data. 


V. SUMMARY 


In this thesis we set out to develop two algorithms for modeling unknown systems 
as ARMA processes. These are the IV method of system identification presented in 
Chapter 3, which is a modification of the method of ordinary least-squares, and the it- 
erative multichannel method presented in Chapter 4. 

The IV method is not a new concept in either its block or sequential processing 
forms. However, our derivation of the sequential algorithm was done independently of 
other IV sequential algorithms. We chose the IV method because it reportedlv has good 
noise handling capabilities and yields consistent and unbiased estimates of the unknown 
svstem’s parameters. It also remains as easy to use as the method of ordinary least- 
squares. We also wanted to gain famiharity with it because it was a possible candidate 
for incorporation into the multichannel method. 

Operating alone, the IV method produces accurate estimates of the unknown sys- 
tem’s parameters. Convergence was within 20 iterations for several second-order sys- 
tems that we tested. Convergence slows down as the system order increases. However, 
the results do converge to the actual system parameters given sufficient processing time. 
The performance of the IV method 1s similar to the performance of the method of ord- 
narv least-squares. 

The proposed iterative multichannel algorithm 1s new in both its block and sequen- 
tial processing forms to the best of our knowledge. It is very simple to use in the block 
form. It achieves accurate results for second-order systems but worse results for third- 
order systems with block correlation elements calculated based on only 500 data points. 
Implementing the stopping parameter increases the accuracy when the parameter esti- 
mates converge but not to the true parameters. Due to its ability to separately process 
the input and output data from the unknown svstem, correlation matrices in the multi- 
channel block processing case are half the size of correlation matrices required for the 
single channel block processing case. This feature reduces the computational intensity 
over what 1s required for the conventional least-squares processing case. The number 
of iterations required for convergence seems to be independent of the order of the svs- 
tem. However, the accuracy of the estimates suffer as the order of the system increases. 
Using linear prediction to estimate the initial values of the AR parameters did not speed 


up convergence or increase the accuracy of the parameter estimates. However, using 
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modified linear prediction significantly increased the accuracy of the AR parameter es- 
timates, although it had no effect on the MA parameter estimates. 

We formulated the multichannel sequential algorithm. This allows the estimates of 
the parameters to be updated as new data becomes available. But we have not tested 
this algorithm. It needs checking using a variety of second and third-order test cases to 
verify that it works. During testing, guidelines need to be developed for the best meth- 
ods to initialize the P and Q matrices to achieve the quickest convergence and the most 
accurate parameter estimates. 

As mentioned above, one of the reasons for investigating the IV method of system 
identification was for possible inclusion into the multichannel algorithm, the hope being 
that the favorable noise performance of the IV method would improve the performance 
of the muluchannel method. This is another area that remains unexplored. 

The block multichannel and IV methods achieved simular results for second-order 
test cases. Convergence to the actual svstem parameters came within 20 iterations for 
both algorithms. However, for third-order systems, convergence was much quicker with 
the multichannel block algorithm than with the sequential IV algorithm. But the pa- 
rameter estimates by the IV method were more accurate than by the multichannel block 
method. A combination of the two algorithms has the potential for incorporating their 
unique advantages into a better overall parameter estimation method. 


Areas for further research are listed below: 


e Verify that the multichannel sequential algorithm developed in Chapter 4 works as 
a means of modeling an unknown system as an ARMA process. 


e Investigate the possibility of incorporating the IV method into the multichannel 
sequential algorithm. 


e Analvze why initializing the AR parameters to the values calculated bv linear pre- 
diction improves the speed of convergence of the AR parameters in the multi- 
channel block algorithm but does not improve the convergence of the MA 
parameters. Identify a method for obtaining an initial estimate of the MA param- 
eters to improve their speed of convergence and accuracy. 


e Investigate the effects of increasing the number of data points used to calculate the 
correlation matrices for the multichannel block algorithm on the accuracy of the 
parameter estimates and their speed of convergence. 


e Investigate the performance of the IV method with noise present. Compare this 
performance with the performance of the method of ordinary least-squares with 
noise present. 
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APPENDIX 


A. INSTRUMENTAL VARIABLE ALGORITHM 


vevtedvededevedescededetenheteKdededesnKkdhkqTeawhKkRccKRkeR RAR «Rca RRR RR a: a: 


> 
we 
er 
ve 


whe 
es 


PAUL DAL SANTO ve 

IV ALGORITHM “ 

4/12/88 ? 

THIS PROGRAM CALCULATES THE AR AND MA PARAMETERS OF A * 
TEST SYSTEM BASED UPON ITS INPUT AND OUTPUT DATA * 
be 


BY USING THE SEQUENTIAL IV METHOD. 


vevededevevevededsedededesvedesedesetokdeleseseledesesosesesesestededesesedesesesestesededtetesedketesteiesetekdesckdesesesedetessetoieve 


aN @¢h eu ee 


VARIABLE DEFINITIONS 


RIVCOF 


QMAT 
IV 

COR 
COEF 
SCALAR 
SCALR2 
SEED 
WSIZE 
YSIZE 
USIZE 


seseses’k 


ARRAY FOR STORING THE AR AND MA PARAMETERS 

THE PROGRAM CALCULATES 

VECTOR OF DATA FROM THE OUTPUT OF THE AUXILIARY 
MODEL AND THE INPUT TO THE TEST SYSTEM 

TRANSPOSE OF VECTOR Z 

VECTOR OF DATA FROM THE OUTPUT AND INPUT OF THE 

TEST SsysiEn 

TRANSPOSE OF THE X VECTOR 

STORAGE FOR INPUT DATA 

STORAGE FOR OUTPUT OF TEST SYSTEM 

STORAGE FOR OUTPUT OF THE AUXILIARY MODEL 

THE Q MATRIX OF THE IV ALGORITHM 

VECTOR OF PARAMETERS CALCULATED BY THE ALGORITHM 
RESULT OF INTERMEDIATE STEP IN ALGORITHM CALCULATION 
VECTOR OF TRUE PARAMETERS OF TEST SYSTEM 

RESULT OF SCALAR INVERSION IN INTERMEDIATE STEP COR 
INTERMEDIATE STEP WHEN CALCULATING THE NEW IV VECTOR 
INITIALIZATION FOR IMSL GAUSSIAN ROUTINE 

ORDER OF AR PART OF THE AUXILIARY MODEL 

OKDEK OF THE AR FART OF THEVIESsles sin 

ORDER OF THE MA PART OF THE TEST SYSTEM AND AUXILIARY 
MODEL 


VARIABLES THAT END IN R CONTAIN THE ROW SIZE OF THEIR RESPECTIVE 


MATRICES. 


VARIABLES THAT END IN C CONTAIN THE COLUMN SIZE OF 


THEIR RESPECTIVE MATRICES. 


VARIABLE DECLARATIONS 


Seveveve 


COMMON /D/ RIVCOF(0: 1000, 10) 


REAL Z(10510),2 TPOC105 Oe xC10 510 he reo erer Ga) 
REAL UG#) 5 _YC10;, 10) Ww GlOr 0) 


4§ 






IVA0015 
IVA0016E 


IVA0025 
IVA0026 
IVA0027 
IVA0028 
IVA0029: 
TVA0030. 
IVAO0031 
IVA0032. 


IVA0044. 
IVA0045° 
IVA0046: 
IVA0047 
IVA0048 


+ id 


Peeeeernt( 10,10) ,0 TEMP( 10, 10) QREMP2( 10,10) ,QTEMP3( 10,10) 
REAL 1V(10,10), 1VTEMP( 10,10) ,COR( 10,10) ,COEF( 10,10) 
REAL SCALAR ,SCALR2 


DOUBLE PRECISION SEED 


INTEGER 1,J,K,WSIZE,YSIZE,USIZE 

INTEGER ZMR,ZMC,ZTR,ZTC,XMR,XMC,XTR, XTC 

INTEGER UMR,UMC, YMR, YMC ,WMR,WMC 

INTEGER QMR,QMC,QTR,QTC,QT2R,QT2C ,QT3R,QT3C 
INTEGER IVR,IVC,IVTR,IVTC,CMR,CMC,COEFMR,COEFMC 
INTEGER IVCR,IVCC 

LOGICAL EOF 


READ(4,*,END=22) YSIZE,USIZE ,COEFMR,COEFMC,ITERA 
CALL RDMAT( COEF , COEFMR , COEFMC) 


INITIALIZE VARIABLES 


EOF = . FALSE. 
ATR = COEFMC 
XTC = COEFMR 
ZTR = COEFMC 
ZTC = COEFMR 
IVR = COEFMR 
IVC = COEFMC 
QMR = COEFMR 
QMC = COEFMR 
SEED = 888. OD1 
IVCR = 0 

IVCC = COEFMR 


WSIZE = YSIZE 


* ZERO OUT THE IV PARAMETER VECTOR, THE AUXILIARY MODEL DATA VECTOR 


AND THE TEST SYSTEM DATA VECTOR. 
CALL INITCIV,IVR,IVC,0. 0) 
CriiielNit¢Z TPO,ZTR,ZIC,0. 0) 
Colieeibtt xX LPO, XTR, XTC, 0: 0) 


INITIALIZE THE QMAT AS A DIAGONAL MATRIX WHOSE DIAGONAL ELEMENTS 


* EQUAL 100 


CALL INITD(QMAT,QMR,QMC, 100. 0) 

GET VALUE FOR U(0). U IS A GAUSSIAN RANDOM VARIABLE. 
CALL GGNML (SEED,1,U) 

SHIFT U(0) INTO X & Z VECTORS TO CREATE X(0) & Z(0) 
VG 0.0 
mies TPO, XTC,YSIZE,USIZE, Y(1,1),,UC1)) 


WO1,1)°>= 0.0 
CALL SHipiez, TPO, ZTCyWSIZE, USIZE WoL) 15 UC1)) 
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IVA00490 
IVA00500 
IVA00510 
IVA00520 
IVA00530 
IVAC0540 
TVA00550 
IVA00560 
IVA00570 
IVA00580 
IVA00590 
TVA00600 
IVA00610 
TVA00670 
TVA00680 
IVA00690 
IVA00700 
IVA00630 
IVA00640 
IVA00650 
IVA00710 
IVA00720 
IVA00739 
IVA00740 
TIVA00750 
IVA00760 
IVA00770 
IVA00780 
IVA00790 
IVA00800 
IVA00810 
IVA00820 
IVA00830 
IVA00840 
IVA00850 
IVA00860 
IVA008/70 
IVA008&0 
IVA00890 
IVA00900 
IVA009 10 
IVA00920 
IVA00930 
IVA00940 
IVA00950 
IVA00960 
IVA00970 
TVA00990 
IVA01000 
IVA01010 
TVA01020 
IVA01030 
IVA01040 
IVA01050 
IVA01060 


* CALCULATE Y(O) = X TPO(O) * COEFFICIENT VECTOR 
CALL MULTICX TPO,XTR,XIC, COEF SCOEFMR , CORP een 1) 
* CALCULATE W(0) = Z TPO(O) * IV VECTOR 


CALL MULTI(Z TPO,ZTR,ZTC,IV,IVR,IVC,W,1,1) 
DO 90 I = O,ITERA 

CALL GGNML(SEED,1,U) 

CALL TPOSE(IV,IVR,IVC,IVTEMP,IVTR, IVTC) 


* SAVE THE IV PARAMETERS 


DO 91 J = 1,IVTC 
RIVCOF(I,J) = IVTEMP(1,J) 
91 CONTINUE 
CALL PRMAT(IVTEMP, IVTR, IVTC) 


* SHIFT Y(M) AND U(M+1) INTO X TPO(M) TO GET X TPO(M+1) 
CALL) SHIFT(X DTPO,XTC,YSIZH.USIZE,Y(1,1),UCt 


* SHIFT WCM) AND U(M+1) INTO 2 TPO(M) TO GET Z TPOG@ TES 
CALL SHIFTI(Z TPO,ZTC,WSIZE,USIZE ,WC1,1), UGS 


aS 


* CALCULATE Y(M+1) AND Z(M+1) 
CALL MULTI(X TPO,XTR,XTC,COEF,COEFMR,COEFMC,Y,1,1) 
CALL MULTIC(Z TPO. ZTREZ TOR) clVRe i cyte 


aS 


* CALCULATE THE NEW Q MATRIX 


CALL MULTI(X TPO,XTR,XTC,QMAT,QMR,QMC,Q TEMP,QTR,QTC) 

CALL TPOSE(Z TPO,ZTR,ZTC,Z,ZMR,ZMC) 

CALL CORE(CMAT,QMR,QMC,Z,ZMR,ZMC,X TPO,XTR,XTC,COR,CMR,CMC) 
CALL MULTI(COR,CMR,CMC,Q TEMP,QTR,QTC,QTEMP2 ,QT2R,QT2C) 
CALL SUBTRC(QMAT,QMR,QMC,QTEMP2 ,QT2R,QT2C, QMAT,QMR,QMC) 


* CALCULATE THE NEW IV VECTOR 


CALL MULTI(X TPO,XTR,XTC,IV,IVR,IVC,IVTEMP, IVTR, IVTC) 
SCALR2 = Y(1,1) - IVTEMP(1,1) 

CALL SMULTI(SCALR2 ,COR,CMR,CMC) 

CALL ADD(IV,IVR,IVC,COR,CMR,CMC,IV,IVR,IVC) 


90 CONTINUE 


* PLOT THE IV PARAMETERS VS THE ITERATION NUMBER 
CALL PLOT2( ITERA, USIZE,YSIZE ,COEF) 


ZZ STOP 
END 


rh sedeteacaededeakecledcakvaesdededesedede deve seskesese tevede sesededcdededeseveve se vesesesese ved desk dese dese dese sete ce seve eek ae 


SUBROUTINE CORE(MAT1,11R,11C,MAT2,12R,12C,MAT3,13R,13C,RMAT, IRR, 
+IRC) 


se tevetevededetedsisdsisditecdk doiedtestehetetietcivietnsekicnkeKildioliicekkckddksddakiheREsted 
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REAL MAT1(10,10) ,MAT2(10, 10) ,MAT3(10,10),RMAT( 10,10) 
REAL Q TEMP(10,10) ,QTEMP2( 10,10) 
INTEGER IRR,IRC,QTR,QTC,QT2R,QT2C 


* CALCULATE THE CORE: Q(M)Z(M+1)°1+X'(M+1)Q(M)Z(M+1) %*-1 
* MAT1 IS THE Q MATRIX, MAT2 IS THE 2 VECTOR, AND MAT3 IS THE 
* X VECTOR. 


CALL MULTI(MAT1,11R,11C,MAT2,12R,12C,Q TEMP,QTR,QTC) 
CALL MULTI(MAT3,13R,13C,Q TEMP,QTR,QTC,QTEMP2 ,QT2R,QT2C) 
SCALAR = 1/(1 + QTEMP2(1,1)) 

CALL SMULTI(SCALAR,Q TEMP ,QTR,QTC) 

CALL ADD(Q TEMP,QTR,QTC,Q TEMP,QTR,QTC,RMAT, IRR, IRC) 
CALL SMULTI(0. 5,RMAT, IRR, IRC) 

RETURN 

END 


x Seicvesoekickcskvedakkkekkkskakdescdvedk ved akedest 
SUBROUTINE PLOT2( ITERA, ICURVN, ICURVD,MAT1) 


Vs Seuicavvevevedescvedecevesesesesleveseveveseveseveseveve veces 


COMMON /D/ RIVCOF(0: 1000, 10) 
REAL X(0: 1000), Y(0: 1000) ,MAT1( 10,10) ,MAX,MIN 
PNG@EGER I,J, LTERA,LCURVN, ICURVD, ISTP 


CALL LIMITSCICURVN, ICURVD ,NMAX,NMIN,NSTP, 
+DMAX , DMIN,DSTP, ITERA) 


* GENERATE THE ITERATION NUMBER 


DO 90 I = 0,ITERA 
X(I) =I 
Y(I) = 0.0 


90 CONTINUE 


© CALCULATE X AXIS LABELING INTERVAL 
ISTP = ITERA/10 


* SET UP THE PLOT 


CALL SHERPA('IVGRAF ','‘A'‘,3) 
CALL RESET(' ALL’ ) 

CALL PAGE(8.5.11.0) 

CALL HEIGHT(0. 14) 

CALL HWROT( ‘AUTO’ ) 

CALL XINTAX 

CALL YAXANG(0) 

CALL PHYSOR(1.5,6. 0) 

GALL AREAZD( Sa0ms. 5) 

CALL COMPLEX 

CALL XNAME('ITERATIONSS' , 100) 
CALL YNAME('COEFFICIENT VALUESS' ,100) 
CALL MESSAG('(A)$' ,100,2.4,-0. 8) 
CALL THKFRM(0. 03) 

CALL FRAME 
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TVA01970 
IVA02000 
IVA02010 
IVA02020 
IVA02030 
IVA02040 
IVA02060 
IVA02070 
IVA02090 
IVA02100 
IVA02110 
IVA02120 
IVA02130 
IVA02140 
IVA02150 
IVA02170 
IVA02180 
IVA02190 
IVA03720 
TVA03730 
IVA03740 
TVA03750 
IVA03760 
IVA03780 
TVA03800 
IVA03810 
IVA03820 
TVA03830 
IVA03850 
IVA03860 
IVA03870 
IVA03880 
IVA03890 
TVA03900 
IVA03910 
TVA03920 
TVA03930 
IVA03950 
IVA03960 
IVA03970 
TVA03980 
IVA04010 
IVA04030 
IVA04040 
IVA04050 
IVA04060 
IVA04070 
IVA04080 
IVA04090 
IVA04100 
IVA04110 
IVA04140 
IVA04150 
IVA04160 
IVA04170 
IVA04180 


CALL GRAF(0, ISTP, ITERA,NMIN,NSTP,NMAX) 
* PLOT THE NUMERATOR VALUES 


DO 93 J = ICURVD + 1,ICURVN + ICURVD 
DO 94 I = 0,ITERA 
Y(I) = RIVCOF(I,J) 
94 CONTINUE 
CALL CURVE(X,Y, ITERA,0) 
93 CONTINUE 


* PLOT DASHED LINES FOR THE COEFS' TRUE VALUES 
CALL DASH 


* PLOT NUMERATOR COEFS' TRUE VALUES 


DO 95 K = ICURVD + 1,ICURVD + ICURVN 
DO 96 J = 0,ITERA 
Y(J) = MATI(K,1) 
96 CONTINUE 
CALL CURVE(X,Y,ITERA,0) 
95 CONTINUE 
CALL ENDGR(0) 


* SET UP SECOND PLOT FOR DENOMINATOR PARAMETERS 


CALL RESET('DASH' ) 
CALIMPHYSO@R@1s 5 01/5) 

CALL SARE Sais. 0,3. 5) 

CALL COMPLX 

CALL XNAME( 'ITERATIONSS' ,100) 

CALL YNAME('PARAMETER VALUESS' ,100) 
CALL MESSAG('(B)S$' ,100,2.4,-0. 8) 

CALL THKFRM(0. 03) 

CALL FRAME 

CALL GRAF(0,ISTP,ITERA,DMIN,DSTP,DMAX) 


* PLOT THE DENOMINATOR VALUES 


DO 91 J = 1,1ICURVD 


DO 92 I = O,ITERA 
YC1)3= -KIVCORGIAA 
oe CONTINUE 


CALL CURVE(X,Y,ITERA,0) 
Z)ih CONTINUE 


% PLOT DENOMINATOR COEFS' TRUE VALUES 


CALL DASH 
DO 99 K = 1,ICURVD 
DO 100 J = 0,ITERA 
Y(J) = -MAT1(K,1) 
100 CONTINUE 
CALL CURVE(X,Y, ITERA,0) 
99 CONTINUE 
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CALL ENDPL(0) 
CALL DONEPL 
RETURN 

END 


Tevedsvessvedessseakdedese sess desea se weve seve seve ss fede Heese sevevevek seve vessscie ak se vesesevesededesesevesesestestesese 
eee ND SSS PTR, LC ieee et 20 , RMAT » IRR, IRC) 


sescsksecestsescsestestesesteseseseses sesrse oat seve Soca rie) ei et eseok> ge ee esisvesesese nn sevevsteveseseveveve 


PURPOSE - ROUTINE SUBTRACTS MAT2 FROM MAT1 AND PUTS THE 
RESULT IN A THIRD MATRIX. 


Reena 1 (10,10), MAT2( 10mg), RMATC 10,10) 
INTEGER IRR, IRC 


Solel sMUie1(<1-0,MAT2.12R.12C) 

CALL ADD(MAT1,11R,11C,MAT2,12R,12C,RMAT,IRR,IRC) 
IRR = I1R 

Ee) SAG, 

RETURN 

END 


ciety se lve le paeh gel oe beasea pela lor tor soy hor oe boar bake hoy ei tea tre ese oe hae ah ae hee sk aoe ee oie ee) aed alae foe oer Use owe ohh oe irae ar ie) eee ok der 3 


ete TNE TPOSE(MAT1, is LIC. ia PRR SRG) 


Sea aes rien) origrigk yrs tsviseyy aris cre lerksr orig eile lat es Fiche isl eee eee. aeveredededevese 
PURPOSE-SUBROUTINE TRANSPOSES A MATRIX AND PUTS THE RESULT 
INTO A NEW MATRIX 


REAL MAT1( 10,10) ,RMAT(10,10) 
Proecek I,J, IRRSIRC 


DORIS I=., 116 

DO 94 J=1,11R 
RMAGTCI , J} 

CONTINUE 

CONTINUE 

IRR = [1c 

IRC = I1R 

be LURN 

END 


= MAT1(J,1) 


MULTICHANNEL ALGORITHM 


sevevevevevetevovededtovedededevedete tote tedectete cede te tote Fede eke ee ERICEIRA REE REERERREABER GS 


ve 

* PAUL DAL SANTO 8/15/88 

sk 

* TWO-CHANNEL SYSTEM IDENTIFICATION ALGORITHM 
¥ 

* PROGRAM CALCULATES THE AR AND MA PARAMETERS 
% BASED UPON THE INPUT AND OUTPUT DATA OF A 

7 TESt- cio LEM, 

ac 

4 SHIFT USED TO CALC RYY FOR THE LINEAR 


oe) 


se 
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EGC) a Ga) 


¥ THE COEFF DATA FILE x 

rd ry 

7% NUMBER OF ITERATIONS IS READ FROM COEFF * 

# DATA FILE * 

oy * 

* we 

a rh 

WRB RT TTR RTT TT Fede TRE RB RITTER ARERR TERERERKKRKKKR IR 

WI VARIABLE DEFINITIONS KK 

RAWDAT ARRAY WHICH CONTAINS THE INPUT AND OUTPUT DATA 
OF THE SYSTEM UNDERSITESE 

Bl ARRAY FOR STORING THE STOPPING PARAMETER AND THE 
COEFFICIENT ERROR 

ARRD ARRAY FOR STORING THE AR PARAMETER ESTIMATES 

ARRN ARRAY FOR STORING THE MA PARAMETER ESTIMATES 

RYYM AUTOCORRELATION MATRIX OF THE OUTPUT DATA 

RYUMARUA CROSSCORRELATION MATRICES 

RUUM AUTOCORRELATION MATRIX OR THE INPUT DATA 

RUUINV INVERSE OF THE AUTOCORRELATION MATRIX OF THE INPUT DATA 

RYYINV INVERSE OF THE AUTOCORRELATION MATRIX OF THE OUTPT DATA 

DM VECTOR OF CURRENT MA PARAMETER ESTIMATES 

CM VECTOR OF CURRENT AR PARAMETER ESTIMATES 

TRUNRM FUNCTION WHICH CALCULATES THE NORM OF THE TRUE VALUES 
OF THE TEST SYSTEM'S PARAMETERS 

WKMAT WORKING MATRIX FOR DOING MATRIX INVERSIONS 

X TPO TRANSPOSE OF THE VECTOR OF INPUT DATA 

af MATRIX FOR THE OUTPUT “OF THE Test eavsote 

U INPUT TO TEE. Teed syste 

COEFM VECTOR OF TRUES@eUEFICIENTS OF THE DEstesSyster 

ITERA NUMBER OF ITERATIONS TO PERFORM 

CORRLN LENGTH OF CORRELATIONS TO USE TO CALCULATE RYY, RUU 
RYU, AND RUY 

Volar ORDER OF THE AR PART OF THE TEST SYSTEM 

Ua bAle ORDER OF THEY MAS PART (OR MIE” TES. Sys ce 

KTIME DHE (CURRENT STE KA TLON 

SHFT SIZE OF THE STARTING LAG FOR LINEAR PREDICTION OF 
THE AR PARAMETEKS 

SEED INITIALIZATION PARAMETER FOR IMSL ROUTINE WHICH 


OF THE AR PARAMETERS IS READ FROM z 


GENERATES RANDOM GAUSSIAN NUMBERS 


INTEGER VARIABLES THAT END WITH R CONTAIN THE ROW SIZE OF 


A PARTICULAR ARRAY. 


INTEGER VARIABLES THAT END WITH C CONTAIN 


THE COLUMN SIZE OF A PARTICULAR ARRAY. 


ciasinciets 


VARIABLE DECLARATIONS Ferere re 


COMMON /D/ RAWDAT(2,0: 2000) ,E1(2,0: 1000), 
+ARRD(0: 1000 ,5),ARRN(0: 1000,5) 


REAL RYYM(5,5),RYUM(5,5),RUUM(5,5),RUYM(5,5),RUUINV(5,5) 
REAL RYYINV(5,5),DM(5,5),CM(5,5),TRUNRM 


INTEGER RYYR,RYYC,RUUR, RUUC , RYUR ,RYUC , RUYR , RUYC 
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25 
* 


* 
* 


"” 


7 


+c 
rid 


INTEGER DR,DC,CR,CC 


Pei 12,12)5x% 1TP0C10,, 10), ¥(2,2) ,UC1),COEFM( 10,10} 
INTEGER WKR,WKC ,XTR,XTC ,COEFR,COEFC , ERR 


INTEGER ITERA,CORRLN,YSIZE,USIZE ,KTIME, SHFT 
DOUBLE PRECISION SEED . 
TEMPORATY MATRICES FOR PERFORMING CALCUALTIONS 
REAL T1M(5,5),T2M(5,5),T3M(5,5) 
i PCH Miies MiGwnaneneC T3136 


BEGIN MAIN PROGRAM 


READ IN THE SIZE OF THE TEST SYSTEM, THE NUMBER OF ITERATIONS 
TO PERFORM AND THE BEGINNING LAG OF LINEAR PREDICTION OF THE 
AR PARAMETERS 
READ(4S.* END=22) YSIZE,USIZE ,COEFR, COEFC,ITERA, SHIT 
READ IN THE TRUE VALUES OF THE COEFFICIENTS OF THE TEST SYSTEM 
CALL RDMAT(COEFM,COEFR , COEFC ) 
INITIALIZE ROW AND COLUMN SIZES OF AS WELL AS OTHER VARIABLES 


CORRLN = 500 


RUUR = USIZE 
RUUC = USIZE 
RYYR = YSIZE + 1 
RYYC = YSIZE + 1 
T3R = YSIZE 

T3C = YSIZE 

DR = USIZE 

DC =1 

CR = YSIZE + 1 
cc =1 

XTR = COEFC 

XTC = COEFR 

SEED = 888. 0D1 
II =0 

¥(1,1) = 0.0 


EaGieao}y = 0.0 
DIVISR = TRUNRM(COEFM,COEFR,YSIZE,USIZE) 


ZERO OUT THE CORRELATION MATRICES, THE MA PARAMETER VECTOR AND 
VECTOR OF INPUT DATA 


CALL INIT(RYYM,RYYR,RYYC,0. 0) 
CALL INIT(RUUM,RUUR,RUUC, 0.0) 
CALL INIT(RYUM,RYYR,RUUC,0. 0) 
CALL INIT(RUYM ,RUUR,RYYC,0. 0) 


D> 
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CALL INIT(DM,DR,DC,O. 0) 
CALL INIT@ aire. mR XTC .0. 0) 


* INITIALIZE THE AR PARAMETER VECTOR 


CALL INITD(CM,CR,CC,1. 0) 
CALL PRMAT(CM,CR,CC) 


* RUN THE FILTER FOR 2000 TIME STEPS TO GENERATE OUTPUT DATA 


+ 


70 


15 


DO 70 KTIME = 0,2000 


GET VALUE FOR U(K). U IS A GAUSSIAN RANDOM VARIABLE. 
CALL GGNML (SEED,1,U) 
SHIFT U(K) INTO X VECTOR TO CREATE X(K) 
CALL SHIFT(X TPO,XTR,XTC,YSIZE,USIZE,Y(1,1),U(1)) 
CALCULATE VALUE OF Y(K) 
CALL MULTI(X TPO,XTR,XTC,COEFM,COEFR,COEFC,Y,1,1) 


SAVE THE INPUT AND OUTPUT DATA 


RAWDAT(1,KTIME) = Y(1,1) 
RAWDAT(2,KTIME) = U(1) 
CONTINUE 


FORMAT (2X,14,2X,F8.5,2X,F8.5) 


* CALCULATE THE CORRELATION MATRICES 


CALL AUTCOR(1,RYYM,RYYR,RYYC ,CORRLN) 
CALL AUTCOR(2,RUUM, RUUR,RUUC ,CORRLN ) 
CALL CRSCOR(1,2,RYUM,RYYR,RUUC ,CORRLN) 
CALL CRSCOR(2,1,RUYM,RUUR,RYYC,CORRLN) 


* INVERT THE AUTOCORRELATION MATRICES OF THE OUTPUT AND INPUT DATA 


CALL LINV2F(RYYM,RYYR,RYYC, RYYINV,0,WKMAT,ERR) 
CALL LINV2F(RUUM , RUUR , RUUC , RUUINV,0,WKMAT,ERR) 


“ MULTIPLY THE INVERSE OF THE AUTOCORRELATION MATRICES BY THEIR 
* RESPECTIVE CROSSCORRELATION MATRICES 


CALL MULTICRUUINV,RUUR, RUUC ,RUYM,RUUR,RYYC,T1M,T1R,T1C) 
CALL MULTICRYYINV,RYYR,RYYC ,RYUM,RYYR,RUUC ,T2M,T2R,T2C) 


* ESTIMATE THE AR PARAMETERS BY LINEAR PREDICTION 


LE CSHEL GE 1) aN 
CALL CORLA4(DR,CM,CR,CC,T3M,CORRLN, SHFT) 
ENDIF 


WRITE (3,26) II,DM(C1,1),DM(2,1) , DM(Gew UNG) biG.) 
WRITE (3,16) II,CMC1,1),CM(2, 1) CMC GG) or > s 
WRITE (3,*) ' 

CALL, SAVE lea. 11] DiHeDRe ba) 
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* CALCULATE THE COEFFICIENT ERROR 


se 


26 


16 


100 


Coie 2  COEFM, COEPR, CM CR, DM DREIT{YSIZE,USIZE ,DIVISR) 
BEGIN THE ITERATIVE PROCEDURE 


DO 100 II = 1,ITERA 
CALCULATE THE VALUES OF THE PARAMETERS AT ITERATION II 


CALL MULTI(T1M,T1R,T1C,CM,CR,CC,DM,DR,DC) 
CALL MULTI(T2M,T2R,T2C,DM,DR,DC,CM,CR,CC) 
CALL SAVE(1,0,11,DM,DR,DC) 
CALL SAVE(2,0,11,CM,CR,CC) 


CALCULATE THE STOPPING PARAMETER 
CALL ERROR(CM,CR,CC,DM,DR,DC,ITI) 


CALCULATE THE COEFFICIENT ERROR 
CALL ERR2(COEFM,COEFR,CM,CR,DM,DR,II,YSIZE,USIZE,DIVISR) 


WRITE (3,26) II,DM(1,1),DM(2,1),DM(3,1) ,DM(4,1),DM(5,1) 
FORMAT (14,1X,'NUM COEF =' ,5(1X,F9. 6 

Wiebe (3 se eee Onn ),CM(2.1),CM(3 i) ,CM(4,1),CMO5, 1) 
FORMAT (14,1X,'DNM COEF =' ,5(1X,F9. 6)) 

WRITE (3,*) ' ERROR = ',E1(1,11),' COEF ERROR = '7/E1(2,11) 
WRITE (3,*) | 


CM(1,1) = 1.0 


CONTINUE 


* END OF THE ITERATIVE PROCEDURE 


7. =— = = 


* PLOT THE PARAMETER VALUES 


Zo 


we 


91 


CALL PLOTICITERA,DR,CR,COEFM) 


STOP 
END 


Vedevevetevedetotodesetodetesetedetetedede dese tetedele secede de FeFede BRR CRIB TREE I: 


FUNCTION TRUNRM (MAT1,M1R,YSZ,USZ) 


sedevesevevesevetededevededetetededetetetetetetetetetek tote FRI ORE REE RT 


RESte ATIC MIR 
INTEGER YSZ,USZ 


VALUE = 0 
DO 91 I = 1,YSZ + USZ 
VALUE = VALUE + (MATI(1,1))**2 
CONTINUE 
TRUNRM = SQRT(VALUE + 1) 
RETURN 
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25 


90 


oi 
gall 


END 


Sevesevedoosdesccdc desesesesesesedescslvdedeslcslecscslescslke ve dk slesedesecslevevescse dese sles ste vies 


SUBROUTINE AUTCORC IFST,MAT1,M1R,M1C,CORRLN) 


Kevescocseseslcucclesesesestcleotocskcscsescsescvesssedesedesestevevedssedessvessseseseskaese 


THIS SUBROUTINE CALCULATES THE AUTOCORRELATION MATRIX USING 
CORRELATIONS OF SIZE CORRLN 


COMMON /D/ RAWDAT(2,0: 2000) ,E1(2,0:1000), 
+ARRD(0: 1000,5),ARRN(0: 1000,5) 

REAL MAT1(M1R,M1C) 

INTEGER I,J,K,CORRLN 


CALL INITCMAT1,M1R,M1C,0.0) 


CALC CORRELATIONS ALONG THE FIRST ROW OF THE MATRIX. THE 
SHIFT = ABSCNUMBER OF THE COLUMN - 1) 
LENGTH OF CORR = D LENGTH + ABS(1 - THE NUMBER OF THE COLUMN) 


ONCE CORR HAVE BEEN CALCULATED FOR THE FIRST ROW, THEY CAN 
BE COPIED INTO OTHER ROWS. HORIZ DISTANCE OF THE PARTICULAR 
ELEMENT FROM THE MAIN DIAGONAL DETERMINES WHICH CORR TO 
COPY. THIS DISTANCE IS GIVEN BY ABS(ROW NUMBER - COLUMN 
NUMBER). 


CORRELATIONS START 200 POINTS FROM THE BEGINNING OF THE 
DATA TO ELIMINATE THE TRANSIENT OF THE TEST SYSTEM. 


DO 90 J = 1,M1C 
ENDPT = 200 + CORRLN - ABS(1 - J) 
DO 93 K = 200,ENDPT 


MAT1(01,J) = MAT1(1,J) + RAWDATCIFST,K-(1-J))*RAWDATCIFST,K) 
CONTINUE 
MATIC 1,73) = MATICI,J)7 (CORRE Bete BoC) 
CONTINUE 


DO 91 I = 2,M1R 


DO 92 J = 1,MI1C 
MAT1(1,J) = MAT1(1,1+ABS(I-J)) 
CONTINUE 
CONTINUE 
RETURN 
END 


Kacdedesededesesede sevesedesedesedesese sede testes seve desesese se lese de ae de seve ve ve ae 


SUBROUTINE COPY(MAT1 ,M1R,M1C ,MAT2 ,M2R,M2C) 
Kickkvdedededededcdesdescvedededkedededededescscsededede dkedekdcdededkicdeseddedededs 


* THIS SUBROUTINE COPIES A MATRIX INTO A SECOND MATRIX. 


REAL MAT1(M1R,M1C) ,MAT2(M2R,M2C) 
INTEGER ST) 


DO 90 I = 1,M1R 
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900 
90 


* 


x 


DO 900 J = 1,M1C 
MAT2(1,J) = MATI(1,J) 

CONTINUE 

CONTINUE 

M2R = MIR 

M2C = MIC 

RETURN 

END 


Tedededcavdedeakatcdedevescadcsdevedeskededcdevedkcseaeakedede vedesdedevcdededesevedkdedeteve sede 


SUBROUTINE CRSCORCIFST,ISND,MAT1,M1R,M1C,CORRLN) 
Iakdcdesdededesksaedesese dese ds te devese de dese se ved dese dededed devedsdevesedeve dese deve 


* THIS SUBROUTINE CALCULATES THE CROSSCORRELATION MATRIX USING 
* CORRELATIONS OF LENGTH CORRLN 


ri 


96 


a7 


NS. 
94 


COMMON 7/7 KRAWEATG2, 0: 2000) ,E1(2,0: 1000), 


TameCo: 1000, 5) ,ARRNCO: 1000,5) 


REAL MAT1(M1R,M1C) 
INTEGER I,J3,K,CORRLN 


Cri INITCMATISMIR,MIC,0. 0) 


FOR CROSS CORRELATION MUST CALC EACH ELEMENT OF THE CORR MATRIX 
SEPARATELY. CORRELATIONS START 200 POINTS FROM THE BEGINNING 
One pirate INATE THE TRANSIENT Ope VEST SISTEM. 


DO 94 I = 1,M1R 
DO 95 J = 1,M1C 
ENDPT = 200 + CORRLN - ABS(I - J) 
IF (J.GE.1) THEN 
DO 96 K = 200,ENDPT 
MAT1(I,J) = MAT1(1,J) + RAWDAT(IFST,K-(I-J)) 


+**RAWDATC ISND,K) 


CONTINUE 
ELSE 
NOT? Ko = 2007E NUP I. 
ATV) — MATIC(T, J) + KRAWDATCIFST,K)* 


+RAWDAT( ISND ,K+(I-J)) 


CONTINUE 
ENDIF 
MAT1(I,J) = MAT1(1,J)/(CORRLN + 1 - ABS(I-J)) 
CONTINUE 
CONTINUE 
RETURN 
END 


sevedledesteslevecesevedeseves: dedededevevedetotededetedecetede dedetedededesedstedestedede Rk Tete 


SUBROUTINE CORLA4(ISIZE,MAT2 ,M2R,M2C ,MAT3,CORRLN, SHIFTT) 


sestevevesesesededesevedevedede dete tetedevetetededtetetede te de de Feed Fett R REAR EAR EE 


* THIS SUBROUTINE CALCULATES THE CORRELATION MATRIX 
* AND THE CORRELATION VECTOR USED FOR LINEAR PREDICTION OF THE 


7% AR PARAMETERS. 


Ii THEN CALCULATES THE INITIAL ESTINATE OF THE 


* AR PARAMETERS AND PASSES THEM BACK TO THE MAIN PROGRAM IN MAT2. 
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DUA0N3820 
DUA03830 
DUA03840 
DUA03850 
DUA03860 
DUA03870 
DUA03880 
DUA03890 
DUA03900 
DUA039 10 
DUA03920 
DUA03950 
DUA03960 
DUA03970 
DUA03980 
DUA03990 
DUA04000 
DUA04010 
DUA04020 
DUA04030 
DUA04040 
DUA04050 
DUA04060 


ae 


=) al 


90 


X* 


» FS) 
PZ 


me 


ope 


COMMON /D/SKAWDAT(Z,0: 2000) Bl 280s ccuse 


+ARRDCO: 1000 ,5) , ARRNCO: 1000,5) 


REAL MAT2(M2K(N2C) |MAIS( M2R=1 Zhe 
REAL T2M(5,1),93M(5,1),T4N( S25 )ewieat So) 
INTEGER ERR, TIR,T1C,T2R,1T2C, TSK fee sisi Zeon BeeeRn 


GENERATE THE FIRST ROW OF THE RYY MATRIX. 
SHIFTS ARE GREATER THAN THE ORDER OF THE 


NUMERATOR. 
M3R = M2R-1 
M3C = M3R 


CALI ONT TOMATS Hoke f3C ,.0..0) 


DO 90 J = 1,M3C 
ENDPT = 200 + CORRLN - SHIFTT 
DO 91 K = 200,ENDPT 
MAT3(1,J) = MAT3(1,J) + RAWDAT( 1,K+SHIFTT)*RAWDAT( 1,K) 
CONTINUE 
MAT3(1,J) = MAT3(1,J)/(CORRLN-SHIFTT+1) 
SHIFTT = SHIFTT + 1 
CONTINUE 


COPY ELEMENTS FROM THE FIRST ROW INTO OTHER LOCATIONS 


De 82 1 = DIGe 
DO 93 J = 1,M3C 
MAT3(1,J) = MAT3(1,1+ABS(I-J)) 
CONTINUE 
CONTINUE 


GENERATE THE RYY VECTOR BY COPYING ELEMENTS OF THE RYY MATRIX 
T2R = M3C 


T2C = 1 
CALL FILL(1,T2R-1,1,1,T2M,T2R,1,2,1,MAT3 ,M3R,M3C) 


GENERATE THE LAST ELEMENT IN THE RYY CORRELATION VECTOR. 


FINELE = 0.0 
ENDPT = 200 + CORRLN - SHIFTT 
DO 94 K = 200,ENDPT 
FINELE = FINELE + RAWDAT(1,K+SHIFTT)*RAWDAT(1,K) 
CONTINUE 
FINELE = FINELE/(CORRLN+1-SHIFTT) 


COPY THE FINAL ELEMENT INTO THE VECTOR 
T2M(T2R,1) = FINELE 


CALCULATE THE INITIAL ESTIMATE OF THE AR PARAMETERS 


CALL LINV2F(MAT3 ,M3R,M3C ,T4M,0,WKMAT, ERR) 
CALL MULTILCT4M,M3R MSC. 12M, 12k 126, Tan ska Gy 


COPY THE AR PARAMETERS INTO THE RETURN ARGUMENT 
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DUA0415 
DUA0416 
DUA0417. 
DUA0418 
DUA0419 
DUA0420 
DUA0421 
DUA04229 
DUA04 23° 
DUA0424 
DUA0425: 
DUA0426 
DUA0427- 
DUA0428) 
DUA0429: 
DUA0430 
DUA04311 
DUA0432: 
DUA0433' 
DUA0434: 
DUA0435! 
DUA04 36! 
DUA0437! 
DUA04 38' 
DUA0439t4 
DUA0440! 
DUA0441i 
DUA04421 
DUA0443! 
DUA0444t 
DUAO445t 
DUA04S46t 
DUA0447! 
DUA0448i 
DUA04491 
DUAO04 501 
DUA045 11 
DUA0452( 
DUA045 31 
DUA0454( 
DUA0455( 
DUA0457( 
DUA0458¢ 
DUA0463( 
DUA0464( 
DUA0465( 
DUA0467( 
DUA0469( 
DUA04 7 0¢ 
DUA047 1( 











ve 


ae 


MATZCi 1) = 1.0 
DO 95 L=1 
Wee2(L + Ls orien. |) 
MAT2(L,1) = T3MCi-1,1) 
Wie (3,*) MAT2ZCL, 1) 
CONTINUE 


Le 
: 
1, 


RETURN 
END 


dedededededededededevedevedededese de deve dedeskedesdkedevedeseseskcsesedesede ses dese ses seve ve 


SUBROUTINE ERROR(MAT1,M1R,M1C ,MAT2 ,M2R,M2C, ITNUM) 
deaesedededededededesevevedevevedededededevedetervedevevesetededededevededeseves fesekeve de 


* THIS SUBROUTINE CALCULATES THE STOPPING PARAMETER. 


ell 


a2 


90 


rg 


COMMON /D/ RAWDAT(2,0: 2000) ,E1(2,0: 1000), 


Teme CO: 1000.5) ,ARRNCUO: 1000 ,5) 


REAL MAT1(M1R,M1C) ,MAT2(M2R,M2C) 
ENCEGER I,J,k 


ERVAL = 0.0 
XVAL = 0.0 
ZVAL = 0.0 


DO 90 I = 400,450 
DO 91 J = M1R,1,-1 
XVAL = XVAL + MAT1(J,1)*RAWDAT( 1, I1+M1R-J) 
CONTINUE 
DO 92 K = M2R,1,-1 
ZVAL = ZVAL + MAT2(K,1)*RAWDAT( 2, I+M2R-K) 
CONTINUE 


ERVAL = ERVAL + (XVAL-ZVAL)**2 
XVAL = 0.0 
ZVAL = 0.0 

CONTINUE 


E1(1,ITNUM) = ERVAL/51 


RETURN 
END 


em! aft - 
Yevevesevevovetetevecedededtetetetcdedetedetedesesedotededededesesededsdetedededededevede de 


SUBROULTINE ERR2(MAT1,MIR,MAT2 ,M2R MATS .M3R,ITNUM, 


+YSZ,USZ,DIVISR) 


Sesevetevedesdedescdedededesdesesedevdesesesdedestededesesevedededke teak dededede deve deve vevese se 


* THIS SUBROUTINE CALCUALTES THE COEFFICIENT ERROR. 


COMMON /D/ RAWDAT(2,0: 2000) ,E1(2,0: 1000), 


+ARRD( 0: 1000,5) ,ARRN( 0: 1000,5) 


REAL MAT1(M1R, 1) ,MAT2(M2R, 1) ,MAT3(M3R, 1) 
INTEGER 1,YSZ,USZ 


ERRVAL = (1 - MAT2(1,1))%**2 
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DUA04720 
DUA04730 
DUA04740 
DUA04750 
DUA04760 
DUA04770 
DUAG+600 
DUA04810 
DUA04820 
DUA04830 
DUA04840 
DUA04850 
DUA04860 
DUA04870 
DUA04880 
DUA04890 
DUA0N4900 
DUA04910 
DUA0N4930 
DUA04940 
DUA04950 
DUA04960 
DUA04970 
DUA04980 
DUA04&990 
DUA05000 
DUA05010 
DUA05020 
DUA05030 
DUA05040 
DUA05050 
DUA05060 
DUA05070 
DUA05080 
DUA0N5090 
DUA05100 
DUA05110 
DUA05120 
DUA05130 
DUA05140 
DUAQ5150 
DUA05 160 
DUA05170 
DUA05190 
DUA05 200 
DUA05210 
DUA05220 
DUA05230 
DUA05 240 
DUA05250 
DUA05260 
DUA05270 
DUA05290 
DUA05300 
DUA05320 
DUA05 330 


DUA0534C 





DO 90 1 =sersZ DUA0535C 
ERRVAL = ERRVAL + (MAT1(1,1) + MAT2(I+1,1))**2 DUA0536C 

90 CONTINUE DUA0537¢ 
DUAON5 38¢ 

DO 92 I = 1,USZ DUA0539¢ 
ERRVAL = ERRVAL + (MATIC ItYSZ,1) = MAT3(1,1))**2 DUA0540¢ 

92 CONT INUE DUA0541(¢ 
DUA0542( 

E1(2,ITNUM) = SQRT(ERRVAL)/DIVISR DUA0S43¢ 
DUA0544( 

RETURN DUA0545(¢ 

END DUA0546( 
DUA0N548( 

¥ Jededededevedededevededededesctedk seve Hactedede sede sedesedevete Kies sete ke seseTeH sedeK TeTede sede Ka DUA0549¢ 
SUBROUTINE FILL(I,J,K,L,MAT1,M1R,M1C,R2,C2,MAT2,M2R,M2C) DUAOS50( 

2c savveddeaged de Breeches keene miseedeleisiedeicieieisesiescieveledededededeseedesesedetes DUAO0551(¢ 
DUA0552( 

* THIS ROUTINE FILLS MAT1 FROM MAT2. POSITIONS DUA0553¢ 
* IN MATi FROM (I,K) TO (J,L) ARE FILLED WITH AN EQUAL NUMBER DUA0554( 
* OF ELEMENTSsPROM MAT2 STARTING AY POSITION (R2,C2)- DUA0555( 
DUA0556( 

REAL MAT1(M1R,M1C) ,MAT2(M2R,M2C) DUAOS57( 
INTEGER ROW,COL,R2,C2,C22 DUA0559( 
DUA0560( 

C22 = C2 DUA0N561( 
DUA0N562( 

DO 90 ROW = I,J DUA05 64( 

DOs Cal = KL DUAO565( 
MAT1CROW,COL) = MAT2(R2,C2) DUA05 66( 

G2.= 02 +71 DUA0567( 

OW CONTINUE DUA0568( 
R2 = R2 + 1 DUA0N5 696 

C2soeCZ2 DUA0570( 

SN, CONTINUE DUA0571{ 
DUA0572( 

RETURN DUA05 736 

END DUA0574t 
DUA05756 

ri Tevevevevedede dese develo cede PeR Ve sede KCTS Te PEPE eT TE ESET TOT CTP CETERA DUAO7OO6 
SUBROUTINE Cheeses one CEMA, CESTE » ITERA ) DUAO0701( 

7 Hie arbirtsrtsc srisrfirtir teria nici bro nectar ers rors ete sets TER rari ier oo Pyaar tartar fae torigrtor tris DUAO702¢ 
DUA0703¢8 

ve ROUTINE CALCULATES THE LIMITS FOR THE GRAPH OF THE STOPPING DUAO7041 
* PARAMETER AND THE COEFFICIENT ERROR. DUAO705t 
* STOPPING PARAMETER LIMITS ARE RETURNED IN EMAX AND ESTP. DUA07 06! 
* COEFFICIENT ERROR LIMITS ARE RETURNED IN CEMAX AND CESTP. DUA0707! 
DUA0708' 

COMMON /D/ RAWDAT(2,0: 2000) ,E1(2,0: 1000), DUA0709 
+ARRD(0: 1000,5) ,ARRN( 0: 1000,5) DUAO7 10! 

REAL EMAX,ESIP DUA0712:! 
INTEGER ITERA DUAO7 13! 

DUAO07 14: 
EMAX = 0.0 DUA0715) 

DO 90 I = 0,ITERA DUA0716 

IF (E1(1,1).GT.EMAX) THEN DUAO718. 
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EMAX = E1(€1,1) 
ENDIF 
90 CONTINUE 


EMAX = 1.25 * EMAX 
ESTP = EMAX/5 
CEMAX = 0.0 


DO 91 I = O,ITERA 
IF (E1(2,1).GT. CEMAX) THEN 
CEMAX = E1(2,1) 
ENDIF 
zal CONTINUE 


CEMAX = 1.25 * CEMAX 
CESTP = CEMAX/5 
RETURN 

END 


qv dasdecdedevedededetetesedcsvavdodcdeakvacdesdesede as ae deve de 
SUBROUTINE PLOTICITERA, ICURVN, ICURVD ,MAT1) 


ac Seicveveseddevesevesevededevedevedssevesevetevede de vede seve 


* THIS ROUTINE GENERATES SEPARATE PLOTS OF THE MA 

* AND AR PARAMETERS. IT THEN REPLOTS THE THESE CURVES ALONG 
* WITH THE STOPPING PARAMETER. FINALLY IT PLOTS THE STOPPING 
* PARAMETER AND THE COEFFICIENT ERROR ON THE SAME GRAPH. 


COMMON /D/ RAWDAT(2,0: 2000) ,E1(2,0: 1000), 
+ARRD(0: 1000,5),ARRN(0: 1000,5) 

REAL X(0: 1000) ,¥(0: 1000) ,MAT1(10,10) 

REAL STP,RNMIN,RNSTEP, RNMAX,RITERA 

INTEGER I,J,11R,11C,VAL, ITERA, ICURV 


CALL LIMITS( ICURVN, ICURVD,DMAX,DMIN,DSTEP, 
+NMAX,NMIN,NSTEP, ITERA) 
CALL LIM2(EMAX,ESTP,CEMAX,CESTP, ITERA) 


* GENERATE THE ITERATION NUMBER 


Domo 1 = 0,1TERA 
X(I) = I 
Y(I) = 0.0 


90 CONTINUE 


*% CALCULATE X AXIS LABELING INTERVAL 
STP = ITERA/10. 


* SET THE DEVICE CALLS FOR ALL OF THE PLOTS 


CALL SHERPA('MCGRAF ','A',3) 
CALL RESET( ‘ALL’ ) 


% SECTION 1 


* THIS SECTION GRAPHS THE NUMERATOR AND DENOMINATOR 
* PARAMETERS ON SEPARATE GRAPHS ON THE SAME PAGE. 
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DUAO07190 
DUAO7200 
DUA07220 


DUAO7 230 
DUA07240 
DUA07280 


DUA07290 
DUA07310 
DUA07320 
DUA07330 
DUAO07350 


DUA07360 
DUA07370 
DUA07380 
DUA07390 
DUA07400 
DUA07670 
DUA07680 
DUA07690 
DUA07700 
DUAO7710 
DUAO07720 
DUA07730 
DUA07740 
DUA07750 
DUAO07760 
DUAO7770 
DUA07780 
DUA07790 
DUA07800 
DUA07810 
DUAO07 820 
DUA07830 
DUAO07 840 
DUA07870 
DUA07880 
DUAO07890 
DUAO7900 
DUA07910 
DUAO07920 
DUA07930 
DUA07940 
DUA07950 
DUA07970 
DUA07980 
DUA07990 
DUA08000 
DUA08030 
DUA08050 
DUA08060 
DUA08070 
DUA08080 
DUA08090 


* SET UP THE PLOT FOR THE NUMERATOR COEFF 


CALL PAGE(8.5,11. 0) 

CALL HEIGHT(0. 14) 

CALL HWROT('AUTO' ) 

CALL XINTAX 

CALL PHYSOR(1.5,6. 0) 

CALL AREA2D(5.0,3.5) 

CALL COMPLX 

CALL YAXANG(0) 

CALL XNAME( 'ITERATIONSS' ,100) 

CALL YNAME(' PARAMETER VALUESS' ,100) 
CALL MESSAG('(A)$',100,2.4,-0.8) 
CALL THKFRM(0. 03) 

CALL FRAME 

CALL GRAF(0. ,STP, ITERA,NMIN,NSTEP, NMAX) 


* PLOT THE NUMERATOR PARAMETERS 


DO 93 J = 1,ICURVN 
DO 94 I = 0,ITERA 
Y(I) = ARRN(I,J) 
94 CONTINUE 
CALL CURVE(X,Y,ITERA,0) 
93. CONTINUE 


* PLOT DASHED LINES FOR THE TRUE VALUE OF THE PARAMETERS 


CALL DASH 
DO 97 K = ICURVD,ICURVD+ICURVN-1 
DO 98 J = 0,ITERA 
Y(J) = MAT1(K,1) 
98 CONTINUE 
CALL CURVE(X,Y,ITERA,0) 
97 CONTINUE 
CALL ENDGR(0) 


* SET UP THE PLOT FOR THE DENOMINATOR PARAMETERS 


CALL RESET('DASH') 

CALEY PHYSOR@ie5 41. 5) 

CALL ARE AZT @Sm0ns 95) 

CALL XNAME( 'ITERATIONSS' , 100) 

CALL YNAME( 'PARAMETER VALUESS' ,100) 
CALL MESSAG('(B)$',100,2.4,-0. 8) 

CALL THKFRM(0. 03) 

CALL FRAME 

CALL GRAF(0. ,STP, ITERA, DMIN, DSTEP, DMAX) 


* PLOT THE DENOMINATOR PARAMETERS 


DO 95 J = 1,ICURVD 
DO 96 I = O,ITERA 
Y(I) = ARRD(I,J) 
96 CONTINUE 


DUA0825C 
DUA0826C 
DUA0827( 
DUA08 28C 
DUA0829C 
DUA0832¢ 
DUA0834C 
DUA0835C 






DUA0844C 
DUA0845C 
DUA0846C 
DUA08470€ 
DUA08480 
DUA08490 
DUA08500 


DUA08530 
DUA08540 
DUA08559 
DUA08560 
DUA08570 
DUA08590 
DUA08600 


DUA08640 
DUA08650 
DUA08660 
DUA08680 
DUA08700 
DUA08710 
DUA08720 
DUA08730 
DUA08740 
DUA08750 
DUA08760 


CALL CURVE(X,Y,ITERA,0) 
95 CONTINUE 


* PLOT DASHED LINES FOR THE TRUE VALUES OF THE DENOM PARAMETERS 


CALL DASH 
DO 99 K = 1,ICURVD-1 
DO 100 J = 0,ITERA 
Y(J) = -MAT1(K,1) 
100 CONTINUE 
CALL CURVE(X,Y, ITERA,0) 
99 CONTINUE 


DO 105 J = 0,ITERA 
Y(J) = 1.0 
105 CONTINUE 
CALL CURVE(X,Y, ITERA,0) 
CALL ENDPL(0) 


SECTION: 2 
* THIS SECTION PUTS THE STOPPING PARAMETER ON THE 
* GRAPHS OF THE NUMERATOR AND DENOMINATOR PARAMETERS. 


* SET UP THE PLOT FOR THE NUMERATOR PARAMETERS 


CALL RESET( 'DASH') 

CALL HWROT(' AUTO’ ) 

CALL XINTAX 

CALL PHYSOR(1.5,6. 0) 

CALL AREA2D(5.0,3.5) 

CALL COMPLX 

CALL YAXANG(0) 

CALL XNAME( 'ITERATIONSS’ ,100) 

CALL YNAME('PARAMETER VALUESS' ,100) 
GAME timScAG( (A)S 7100,2.4.-0. 8) 
CALL THKFRM(0. 03) 

CALL FRAME 

CALL GRAF(0. ,STP,ITERA,NMIN,NSTEP , NMAX) 


* PLOT THE NUMERATOR PARAMETERS 


DO 200 J = 1,1CURVN 
DO 201 I = 0,ITERA 
Y(I) = ARRN(I,J) 
Oat CONTINUE 
CALL CURVE(X,Y,ITERA,0) 
200 CONTINUE 


* PLOT DASHED LINES FOR THE TRUE VALUES OF THE PARAMETERS 


CALL DASH 
DO 202 K = ICURVD, ICURVD+ICURVN-1 
DO 203 J = 0,ITERA 
Y(J) = MAT1(K,1) 
203 CONTINUE 


DUA08770 
DUA08 780 
DUA08790 
DUA0D8800 
DUA08810 
DUA08&820 
DUA08830 
DUA08840 
DUA08850 
DUA08 860 
DUA08870 
DUA08880 
DUA08890 
DUA0N8900 
DUA08910 
DUA08920 
DUA08930 
DUA08940 
DUA0 8950 
DUA08960 
DUA08970 
DUA08980 
DUA08990 
DUA09010 
DUA09020 
DUA09030 
DUA09040 
DUA09050 
DUA09060 
DUA09090 
DUA09100 
DUA09130 
DUA09140 
DUA09150 
DUA09160 
DUA09170 
DUA09180 
DUA09190 
DUA09210 
DUA09220 
DUA09230 
DUA09240 
DUA09250 
DUA09260 
DUA09270 
DUA09280 
DUA09290 
DUA09300 
DUA09310 
DUA09320 
DUA09330 
DUA09340 
DUA09350 
DUA09 360 
DUA09370 


CALL CURVE(X,Y,ITERA,0) DUA09380 


202 CONTINUE DUA09390 
DUA09400 

DUA09410 

* PLOT THE STOPPING PARAMETER ON THE SAME GRAPH DUA09420, 
DUA09430 

CALL DOT DUA09440 

CALL YGRAXS(0. 0,ESTP,EMAX,3.5,'STOPPING PARAMETERS’, DUA09450 
+-100,5.0,0.90) DUA09460 
DUA09470 

DO 204 J = 0,ITERA DUA09480 

Y(J) = E1(1,J) DUA09490 

204 CONTINUE DUA09500 
CALL CURVE(X,Y,ITERA,0) DUA09510 

CALL ENDGR(0) DUA09520 
DUA09530 

DUA09540 





* SET UP THE PLOT FOR THE DENOMINATOR PARAMETERS DUA09550 
DUA09560 

CALL RESET( 'DOT’) DUA09570 

CAI ee SO @ieeseleeS ) DUA09580 

CALL AREA2D(5.0,3.5) DUA09590 

CALL XNAME('ITERATIONSS' ,100) DUA09610 

CALL YNAME(' PARAMETER VALUESS' ,100) DUA09620 

CALL MESSAG( '(B)S',100,2.4,-0. 8) DUA09630 

CALL THKFRM(0. 03) DUA09640 

CALL FRAME DUA09650 

CALL GRAF(0. ,STP,ITERA,DMIN,DSTEP, DMAX) DUA09660 
DUA09680 

* PLOT THE DENOMINATOR PARAMETERS DUA09690 
DUA09700 

DO 2055 — 91, FOUR VD DUA09710 

DO 206 I = 0,ITERA DUA097 20 

wel) = ARRDCIaen DUA09730 

206 CONTINUE DUA09740 
CAL IMeURVE( X,Y, I TERA OD DUA09750 

205 CONTINUE DUA09760 
DUA09770 

* PLOT DASHED LINES FOR THE TRUE VALUES OF THE PARAMETERS DUA09780 
DUA09790 

CALL DASH DUA09800 

DO 207 K = 1,ICURVD-1 DUA09810 

DO 20sec) — soe MERA DUA09820 

Y(J) = -MAT1(K,1) DUA09830 

208 CONTINUE DUA09840 
CALL CURVE@ aie TERA.) DUA09850 

207 CONTINUE DUA09 860 
DUA09870 

DO 209 J = 0,ITERA DUA09880 

Y( J) =si0 DUA09890 

209 CONTINUE DUA09900 
CALL CURVE(X,Y,ITERA,0) DUA09910 
DUA09920 


* PLOT THE STOPPING PARAMETER ON THE SAME GRAPH DUA09930 
DUA09940 
CALL DOT DUA09950 
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CALL YGRAXS(0. 0,ESTP ,EMAX,3.5, ‘STOPPING PARAMETERS’ , 
tH. 0,0. 0) 


DO 210 J = 0,ITERA 
Y(J) = El(1,J) 
210 CONTINUE 
CALL CURVE(X,Y,ITERA,0) 
CALL ENDPL(0) 


aOECT LONS3 
* THIS SECTION PLOTS THE STOPPING PARAMETER AND THE COEFFICIENT 
* ERROR ON THE SAME GRAPH. 


Secor lUr THESPLOT FOR THE STOPPING PARAMETER 


CALL DOT 

CALL HWROT('AUTO' ) 

CALL PHYSOR(1.5,6. 0) 

CAE ARH eeD( 5.0325) 

CALL XNAME( 'ITERATIONSS' ,100) 

CALL YNAME('STOPPING PARAMETERS’ , 100) 
CALL THKFRM(0. 03) 

CALL FRAME 

CALL GRAF(0. ,STP,ITERA,0. ,ESTP ,EMAX) 


DO 306 J = 0,ITERA 
Y(J) = E1(1,J) 
306 CONTINUE 
CALL CURVE(X,Y,ITERA,0) 


* PLOT THE COEFFICIENT ERROR ON THE SAME GRAPH 


CALL CHNDOT 
CALL YGRAXS(0. 0,CESTP,CEMAX,3.5, COEFFICIENT ERRORS’, 
oO... 0.0 ) 


DO 307 J = 0,ITERA 
Y(J) = E1(2,3) 
307 CONTINUE 
CALL CURVE(X,Y,ITERA,0) 
CALL ENDPL(0) 
CALL DONEPL 
RETURN 
END 


ve Sevededesevededededesedetesesedekck seth tetelekdeictiolnsiok 


SUBROUTINE SAVE(VAL,M,K,MAT1,M1R,M1C) 


¢ re) - a 7 = = =a = 
sc Yededededevedesedetedededs -selekeveisiclietekscietckiolelentdteiok 


* ROUTINE SAVES PARAMETER ESTIMATES IN EITHER ARRD OR ARRN 
* DEPENDING UPON THE VALUE OF VAL. 


COMMON /D/ RAWDATC2,0: 2000) ,E1(2,0: 1000), 
+ARRD( 0: 1000,5) ,ARRNCO: 1000,5) 
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DUA09960 
DUA09970 
DUA09980 
DUA0N9990 
DUA10000 
DUA19010 
DUA10020 
DUA15030 
DUA10040 
DUA10050 
DUA10060 
DUA10070 
DUA10080 
DUA10090 
DUA10100 
DUA10110 
DUA10120 
DUA10130 
DUA10140 
DUA10150 
DUA10160 
DUA10180 
DUA10190 
DUA10200 
DUA10220 
DUA10230 
DUA10240 
DUA10250 
DUA10260 
DUA10270 
DUA10280 
DUA10290 
DUA10300 
DUA10310 
DUA10320 
DUA10330 
DUA10340 
DUA1i0350 
DUA10360 
DUA10370 
DUA10380 
DUA10390 
DUA10400 
DUA10410 
DUA10420 
DUA10430 
DUA10440 
DUA10860 
DUA10870 
DUA10880 
DUA10890 
DUA10900 
DUA10910 
DUA10920 
DUA10930 
DUA10940 


ELSEIF (VAL. EQ. 2) THEN DUA1103¢ 

ARRD(K,I) = MAT1(I,J) DUA1104¢ 

ENDIF DUA1105¢ 

900 CONTINUE DUA1106( 
90 CONTINUE DUA1107¢ 
DUA1108(¢ 

RETURN DUA1109¢ 


REAL MAT1(M1R,M1C) DUA1096( 
INTEGER I,J,K,M,VAL DUA1097¢ 

DUA1098¢ 
DO 90 I = 1,M1R DUA1099¢ 


DO 90072 —s i Mic 
TPC Ales!) THEN 
ARRN(K,I) = MATI1(I,J) 


END 





DUA1102¢ 


C. SUBPROGRAMS COMMON TO BOTH SYSTEM IDENTIFICATION 
ALGORITHMS SUB80010 








SUBO00Z0 
ste sees’ eve s'e Je esleves PrP es23 eve evesle eves yee le teules eaten’ ee Je Je estes eveslestse seve ese sevesk Jeslevesle sess ese Je ses'e Je Jes! SUBOOO4GC 
SUBROUTINE ay (MAT 1, ER? Tel MARZ, IC2 pitas IRR 5 Je) SUBOOOSC 
s'e ses s'e Jove id oid ve se esle eviews esieste evlevevevdes eve 6 Je levecios rs voneuee ates eo slewle <3 Je id Je Pee este Je hd SUBOQO06C 
SUBOOO7C 
¥ THIS SUBROUTINE ADDS TWO EQUAL SIZE MATRICIES AND PUTS THE RESULT SUBOOO8C 
¥ IN A THIRD MATRIX. SUBOOO9C 
SUBOO10C 
REAL. MATIC TR? ,.1C]) HATZ OUR2 C2 >, RAT ORR aG el Ca) SUBO013C 
INTEGER I,J,IRR,IRC SUB0014C 
SUBO016C 
DO 92 a= ee SUB0N01/7C€ 
DO 920 J=1,ICl SUBO018C 
RMAT(I,J) = MAT1(I,J) + MAT2(I,J) SUB0019C 
220 CONTINUE SUBOO020C 
ss CONTINUE SUBOO21C 
IRR = IR1 SUB00220) 
IRC = ICl SUB0023C 
RETURN SUB0024C 
END SUB00250 
SUB00260 
se Je elo ewe Je we alee! fen s'e eslewe Servevesesics levees eve Je le evesles Pers lo ewe aloute a’ evesesleveste * SUB00280. 
SUBROUTINE INETCMATI MIR, MIC, ENITVL) SUB00290. 
as Jecevedede desc deve dece Fess se vee Te TGC PETC Tee TER OTR TIEN SUB00300 
SUB00310 
* THIS SUBOUTINE INITIALIZES A MATRIX TO INITVL SUB00320 
SUB0N0330 
REAL MAT1(M1R,M1C) , INITVL SUB00340 
INTEGER I,J SUB00350 
SUB00360 
DO 94 I=1,M1R SUB00380 
DO 95 J=1,M1C SUB00390 
MAT1(1,J)=INITVL SUB00400 
95 CONTINUE SUB00410 
94 CONTINUE SUB00420 
RETURN SUB00430 
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“al 
90 


END 


Seveskekocsesesvacseskdcaeslese dese sak steaks Jedesesededesesekdeveveve te 
SUBROUTINE INITDCMAT1,M1R,M1C,INITVL) 


svocdesedeseacskvocdedealsscsesescsevesevdesedtoesescocscsevesevesedesde dese 


THIS SUBROUTINE INITIALIZES A MATRIX TO AS A DIAGONAL MATRIX 
WHOSE DIAGONAL ELEMENTS EQUAL INITVL. °* 


REAL MAT1(M1R,M1C), INITVL 
INTEGER I,J 


DO 94 I=1,M1R 
DO 95 J=1,M1C 
IF (I1.EQ. J) THEN 
MAT1(1,J)=INITVL 
ELSE 
MAT1(1,J)=0. 0 
ENDIF 
CONTINUE 
CONTINUE 
RETURN 
END 


Fesesededelesdevedcdesdedededesesdke dededesededesedededesdededcdevesedesededesesde te dedese deve dese deve sex 
SUBROUTINE LIMITS(NSZ,DSZ ,DMAX ,DMIN, DSTEP ,NMAX,NMIN,NSTEP, ITERA) 


Seresesevedeseseseseadedesesedede sede sesetedesesese ceded se reves devede devetesedevevese ses te 46 Fe 
ROUTINE CALCULATES THE LIMITS FOR THE GRAPHS 
CALCULATES DENOMINATOR AND NUMERATOR LIMITS SEPARATELY 
IN PREPARATION FOR MAKING TWO GRAPHS 


COMMON /D/ RAWDAT(2,0: 2000) ,E1(2,0: 1000), 


FARROCO: 1000 ,5),ARRNCO: 1000,5) 


REAL DMAX,DMIN,DSTEP,NMAX,NMIN,NSTEP 
INTEGER DSZ,NSZ 


CALCULATE THE DENOMINATOR LIMITS 


DMAX 
DMIN 


1eg0 
0.0 
DO 90 I = 1,DS8Z 
DO 91 J = 0,ITERA 
IF ((ARRD(J,1)).GT.DMAX) THEN 
DMAX = ARRD(J,1) 
ENDIF 
IF ((ARRD(J,1)).LT.DMIN) THEN 
DMIN = ARRD(J,1) 
ENDIF 
CONTINUE 
CONTINUE 


IF (DMAX.GT.0) THEN 


DMAX = 1.25 * DMAX 
ELSE 
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SUB00440 
SUBO0450 
SUB00470 
SUBO00480 
SUBO0490 
SUBO0500 
SUB00510 
SUB00520 
SUB00530 
SUB00540 
SUBO00550 
SUB00570 
SUB00580 
SUB00590 
SUB00600 
SUB00610 
SUB00620 
SUB00630 
SUB00640 
SUBO0650 
SUB00660 
SUB00670 
SUB00680 
SUB00690 
SUBO0710 
SUB00720 
SUBO00730 
SUB00740 
SUB00750 
SUB00760 
SUBO0770 
SUB00780 
SUB00790 
SUBO0800 
SUBO0820 
SUB00830 
SUB00840 
SUBO00850 
SUBO0860 
SUB00870 
SUBO00880 
SUBO0890 
SUBO0900 
SUBO00910 
SUB00920 
SUB00930 
SUB00940 
SUBO0950 
SUB00960 
SUB00970 
SUBO0980 
SUB00990 
SUBO01000 
SUB01010 
SUB01020 
SUB01030 


38: 
o2 


oy 


* 


sy 


DMAX = 0.0 


ENDIF 
IF (DMIN. GT. 0) THEN 
DMIN = 0.0 
ELSE 
DMIN = 1.25 * DMIN 
ENDIF 


DSTEP = (DMAX - DMIN)/5 
CALCULATE THE NUMERATOR LIMITS 


NMAX = 0.0 
NMIN = 0.0 
DO 92 I = 1,NSZ 
DO 93 J = 0,ITERA 
IF (ARRN(J,1I).GT.NMAX) THEN 
NMAX = ARRN(J,1) 
ENDIF 


IF CARRN(CJ 1). LT. NMIN) aren 
NMIN = ARRN(J,1I) 
ENDE 
CONTINUE 
CONTINUE 


IF (NMAX. GT. 0) THEN 
NMAX = 1.25 * NMAX 


ELSE 
NMAX = 0.0 
ENDIF 
IF (NMIN.GT.0) THEN 
NMIN = 0.0 
ELSE 
NMIN = 1.25 * NMIN 
ENDIF 
NSTEP = ABS(NMAX - NMIN)/5 
RETURN 
END 


wiivkkkdkvkkeckkeekeeikickkkkkikddsasdsscakddsdedkidiacddede desde desea ve 


SUBROUTINE MULTI (MAT1,M1R,M1C ,MAT2 ,M2R,M2C,RMAT,M3R,M3C) 
Tevedededesede sede sted He ede te dete de de dedede Fede RK TK RK IK 


ROUTINE MULTIPLIES TWO MATRICES AND PUT THE RESULT IN A 
THIRD MATRIX. 


REAL MAT1(M1R,M1C) ,MAT2(M2R,M2C) ,RMAT(M1R,M2C) 
INTEGER I,J,K,IRR,IRC 


CALL INITCRMAT,M1R,M2C,0. 0) 


DOR Ae I=1 MIR 
DO 910 J=1,M2C 
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SUB01210 
SUBO01220 
SUB01230 


5UB01370 
SUB0O1386 
SUB01390 
SUB01400 
SUB01410 
SUB01430 
SUB0O1440 
SUB01450 
SUB01460 


SUB01610 
SUB01620 
SUB01630 
SUB01640 


9100 
710 


OZ 
a2 


C301 
o2 


is 


ois 


DO 9100 K=1,M1C 
RMAT(1,J)=RMAT(1,J) +MAT1(1,K)*MAT2(K,J) 
CONTINUE 
CONTINUE 
CONTINUE 
M3R = MIR 
M3C = M2C 
RETURN 
END 


2197) Se LSTA STATIS LIE Le Ses ok eae hak hea sr br ho hori sol ea srive ieee Lari shy 
SUBROUTINE PRMATCMAT1;I1R,11C) 


Seriedevercsesesevevevesesenevesesesescvesescstsosescvesevesk 


SUBROUTINE PRINTS A MATRIX OUT TO THE FILE DEFINED 
AS UNIT 3 


REAL MAT1(10,10) 
INTEGER I,J 


DO 92 I = 1,11R 
VRIMECGoecO2 mONr1C1.J).J = ieilep 
FORMAT (7(2X,F8.5)) 
CONTINUE 
RETURN 
END 


2 scree: pave hae br (orl se (ae licase Lael oe {seis hor (or sei week Oe se oe ise ie 
SUBROUTINE RDMATCMAT1,M1R,M1C) 


ricriseis ake Seri bistlatlprigkterlscd acter tal ar Lok ek bacdok isk oe bk ok 
ROUTINE READS A MATRIX FROM FILE SPECIFIED AS UNIT 4. 


REAL MAT1(M1R,M1C) 
OWE CER wT 5.J 


READ IN MATRIX 


DO 92 I = 1,MI1R 
READ (4,*) (MAT1(1,J),J=1,M1C) 
FORMAT( 10F3. 1) 
CONTINUE 
RETURN 
END 


Sevedevevevesesosedelesesesesesedeieiedek Retohiiedeklehscdoedekhcselededstetde eke: 


SUBROUTINE SHIFT(MAT1,M1R,M1C,DSIZE,NSIZE, OUTDAT, INDAT) 


CRT a erat sth sin ri tee ote tee iit se trish larieciscirtariictir tr ter tor er Lice ie 

ROUTINE SHIFTS NEW INPUT AND OUTPUT VALUES INTO DATA VECTOR 
OF THE TEST SYSTEM. THE OLDEST VALUES ARE LOST TO MAKE 
ROOM FOR THE NEW VALUES. 


REAL MAT1(M1R,M1C) ,OUTDAT( 1), INDAT(1) 
INTEGER J,DSIZE,NSIZE,NSTART 
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SUB01650 
SUB01660 
SUB01670 
SUB01680 
SUB01690 
SUB01700 
SUBO17 10 
SUB01720 
SUB01730 
SUB01740 
IVA04910 
IVA04920 
1VA04930 


IVA04970 
IVA04980 
IVA04990 
IVA05010 
IVA05020 
TVA05030 
IVA05040 
IVA05050 
IVA05060 
TVA05070 
SUB01760 
SUBO1770 
SUB01780 
SUBO1790 
SUB01800 
SUBO1810 
SUB01840 
SUBO1850 
SUB01860 
SUB01870 
SUBO01880 
SUB01890 
SUBO1900 
SUB01910 
SUB0N1920 
SUB01930 
SUB01940 
SUBO1950 
SUB01970 
SUB01980 
SUB01990 
SUBO2000 
SUBO2010 
SUB02020 
SUB02030 
SUBO2060 
SUB02070 


SUB02080 


SUB02090 
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230 
93 


NSTART = DSIZE + 1 

DO892 J = Dehee a1 
MAT1(1,J) = MAT1(1,J-1) 

CONTINUE 

MAT1(1,1) = OUTDAT(1) 


DO 93 J = MIC,NSTART+1,-1 
MAT1(1,J) = MAT1(1,J-1) 

CONTINUE 

MAT1(1,NSTART) = INDAT( 1) 

WRITE(12,19) (MAT1(1,J),J=1,M1C) 

FORMAT (7(1X,F8. 4)) 


RETURN 
END 


Kosseskeseskueveskdeveskveskedese ve desese de sedesleve dete aka ves te te dede seve dese seve 
SUBROUTINE SMULTE (CONST MAIR MIC) 


tesesieslcsesicdeslevevevescdesvvedivesesessdesevedesevovedevededecdk dese ve de ak ve sk ve ak ve 
ROUTINE MULTI. LIES A VARY ee ooo CALAR. 


REAL MAT1(M1IR,M1C) ,CONST 
INTEGER I,J 


DO 93 I=1,MIR 
DO) 9soed ar tale 
MAT1(1,J) = MAT1(1I,J) * CONST 
CONTINUE 
CONTINUE 


RETURN 
END 


SUB02480 


fe 


SUB02100 
SUB02110 
SUB02120 
SUB02130 
SUB02140 
SUB02150 
SUB02160 
SUB02170. 
SUB02180 
SUB02190 
SUB02200 
SUB02210 
SUB02220 
SUB02230 
SUB02240 
SUB02260 
SUB02270 
SUB02280 
SUB02290 
SUB02300 
SUB02310 
SUB02320. 
SUB02350 
SUB02370 
SUB02380 
SUB02390 
SUB02400 | 
SUB02410 
SUB02420 
SUB02430 
SUB02440 
SUB02450 
SUB02460 
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